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ABSTRACT 


The  advent  of  small,  portable  super  microcomputers  has  enabled  the  rapid  analysis  of 
mechanical  signals.  In  addition  to  the  microcomputer's  reduction  in  cost,  their  capabilities 
have  expanded  to  allow  real  time  analysis  of  the  frequency  characteristics  of  vibrating 
structures.  Data  acquisition  modules  allow  the  addition  of  powerful  components  to  a 
microcomputer.  Modules  consisting  of  the  necessary  Analog  to  Digital  converters.  Digital 
to  Analog  converters,  and  the  hardware  necessary  to  support  sixteen  channels  of  input  data 
can  easily  fit  inside  desk  top  microcomputers.  Fast  32  bit  processors,  small  high  volume 
hard  disk  drives  and  large  RAM  capacities  complete  the  package,  enabling  a  microcomputer 
to  perform  complex  dynamic  signal  analysis.  One  application  of  the  super  microcomputer 
is  in  the  field  of  underwater  explosion  shock  testing.  The  combination  of  the  acquisition 
and  analysis  capabilities  in  the  same  device  eliminates  the  requirement  to  record  the  data  to 
an  intermediate  device  prior  to  analysis.  Thus  the  data  can  be  analyzed  within  minutes  of 
the  completion  of  an  experiment,  giving  immediate  results  at  the  testing  site  instead  of  an 
analysis  facility.  Another  application  of  the  data  acquisition  equipped  microcomputer  is  in 
the  analysis  of  the  frequency  response  characteristics  of  structures.  Programming  the 
microcomputer  to  acquire  the  data  from  a  vibrating  structure  in  order  to  determine  the 
dynamic  signature  is  easily  done.  The  Hilbert  transform,  identified  previously  as  a  means 
of  detecting  non-linearities  in  the  frequency  response,  can  also  be  included  in  the  data 
analysis  portion  of  the  same  program.  Thus  a  complete  package  of  signal  analysis 
routines,  including  frequency  response,  power  spectra,  and  system  linearity,  can  be 
programmed  into  the  microcomputer.  Discussed  in  this  study  is  the  method  utilized  to 
validate  a  data  acquisition  equipped  microcomputer  for  such  applications. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  super  microcomputer  is  a  tool  that  enables  the  engineer  or  scientist  to  rapidly 
digest  large  quantities  of  information.  Recent  advances  in  computer  technology  have 
brought  about  a  startling  increase  in  the  computing  power  of  small  "desktop"  systems.  The 
speed  of  the  processor  has  significandy  increased  while  the  cost  have  dropped  dramatically. 
Coupled  with  the  ability  to  store  large  quantities  of  data  in  internal  memory,  the  speed  of 
computers  have  taken  quantum  leaps  forward.  Five  years  ago  random  access  memory 
banks  of  600  kilobytes  (600K)  were  common  for  a  desktop  scientific  computing  system. 
Today  systems  with  10  Megabytes  are  readily  available.  Processor  word  size  has  increased 
from  8  and  16  bit  word  length  to  32  bit  in  the  same  period.  Data  storage  and  retrieval 
systems  have  also  become  smaller  and  less  expensive.  Thus  the  ability  to  process 
information  rapidly  and  store  it  has  made  micro  and  mini  computer  systems  a  valuable  aid 
to  the  engineer. 

One  of  the  uses  of  the  super  microcomputer  is  the  area  of  data  acquisition. 
Experiments  can  be  rapidly  set  up,  the  data  can  be  collected  and  analyzed  in  near  real  time, 
and  the  results  can  be  displayed  in  fraction  of  the  time  than  was  previously  required. 
Evidence  of  the  increase  in  computing  power  is  evident  in  the  increased  ability  of  Analog 
To  Digital  conversion.  Analog  signals,  frequently  voltages  from  a  set  of  remote  sensing 
devices,  can  be  automatically  read  and  recorded.  In  1973  the  sampling  rate  of  20,000 
samples  per  second  was  considered  good  for  a  microcomputing  system  [Ref.  1].  In 
1986  the  super  microcomputer  is  capable  of  sampling  1,000,000  samples  per 
second  [Ref.  2],  The  only  limitation  is  the  throughput  capability  of  the  data  bus.  This  is  a 
50  fold  increase  in  sampling  rate  and  coupled  with  simultaneous  sampling  cards  and 
multichannel  A/D  converters,  super  microcomputers  are  now  able  to  analyze  signals  in 
wide  band  vibration  experiments  and  shock  analysis  applications. 

In  the  field  of  underwater  explosion  testing  the  engineer  has  been  forced  to  record 
shock  data  to  an  intermediate  device  prior  to  analysis.  This  device,  usually  a  wide  band 
tape  recorder,  is  a  costly  addition  to  the  test  equipment.  If  there  was  a  small,  portable 
device  which  could  record  the  data  and  perform  the  analysis  at  the  test  site  the  engineer 
would  have  immediate  results  available.  Not  only  would  this  save  time,  but  it  would 


eliminate  the  need  for  the  intermediate  recording  device.  The  properly  equipped 
microcomputer  would  allow  the  engineer  to  gather  the  information,  display  the  results,  and 
record  them  for  further  analysis.  This  process  would  take  seconds,  where  as  the  present 
alternative  would  take  much  longer. 

In  the  study  of  the  mechanical  signature  of  structures,  there  has  been  a  recent  emphasis 
toward  use  of  the  Hilbert  transform  in  order  to  determine  if  the  frequency  response  of  a 
structure  is  behaving  in  a  non-linear  manner.  There  are  no  dynamic  signal  analyzers 
available  that  include  the  capability  of  performing  the  Hilbert  transform  of  the  frequency 
response.  In  order  to  perform  this  operation  the  data  would  have  to  be  transferred  to  a 
computer  for  further  analysis.  This  is  a  costly  and  time  consuming  process  If  a 
microcomputer  with  a  data  acquisition  package  could  perform  the  dynamic  signal  analysis 
of  mechanical  signatures  from  a  structure,  then  it  could  also  perform  the  Hilbert  transform 
of  the  frequency  response  of  the  structure.  Again,  this  would  result  in  a  significant 
reduction  in  the  time  required  for  data  analysis  and  include  processes  that  are  not  available 
on  the  present  generation  of  dynamic  signal  analyzers. 

B.  PURPOSE 

The  purpose  of  this  study  is  to: 

1 .  Program  the  properly  equipped  microcomputer  to  perform  data  acquisition,  reduction, 
and  analysis  of  underwater  explosion  shock  data. 

2.  Perform  dynamic  signal  analysis  on  mechanical  signatures  from  structures,  using 
both  impact  and  steady  state  excitations. 

3.  Include  in  the  dynamic  signal  analysis  portion  of  the  data  analysis  routines  the  Hilbert 
transform.  This  would  place  a  complete  signal  analysis  package  in  one  device. 
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A.  DIGITAL  SIGNAL  PROCESSING  FUNDAMENTALS 

The  Fourier  transform  is  a  method  of  calculating  the  frequency  spectrum  of  a  time 

signal.  Fourier  stated  that  any  periodic  time  signal  of  fundamental  frequency  J-  could  be 

*o 

transformed  into  the  summation  of  a  series  of  sinusoidal  signals  of  varying  amplitude  and 
frequency.  The  Fourier  series  of  a  periodic  time  signal  of  period  T0  is  then: 


The  Fourier  transform  of  a  signal  is  a  method  of  computing  the  frequency  spectrum 
of  a  non-periodic  infinite  time  signal,  provided  it  satisfies  Dirichlet's  Conditions.  It  is 
defined  as: 

oo 

X(f)  =  T{x(t)}  =  J  x(t)  e  dt  .  [Eqn.  2] 

-oo 

The  finite  time  signal  is  infinite  in  the  frequency  domain,  the  infinite  time  signal  has  a 
finite  bandwidth  of  frequency  components.  The  fast  Fourier  transform  (FFT)  is  a  recent 
mathematical  technique  allowing  the  rapid  computation  of  the  frequency  components  of  a 
time  signal.  Programed  into  a  microcomputer  with  the  ability  to  convert  a  time  signal  into  a 
set  of  digital  samples,  the  FFT  is  a  powerful  tool  that  allows  fast  computations  of  a  systems 
transfer  function. 

The  process  of  digitizing  an  analog  signal  is  to  multiply  it  by  a  pulse  train  at  a  set 
sampling  frequency.  Given  an  analog  signal  of  bandwidth  FmiU,  the  Nyquist  sampling 
frequency  states  that  the  sampling  rate  of  the  analog  signal  must  be  greater  than  or  equal  to 
2  Fmajl  in  order  to  avoid  aliasing  errors,  thus: 

Fs2  2Fmax.  [Eqn.  3] 


10 


Multiplication  in  the  time  domain  is  equivalent  to  convolution  in  the  frequency 
domain.  This  simply  takes  the  frequency  components  of  the  signal  and  repeats  them  in  the 
frequency  domain  at  intervals  of  the  sampling  frequency.  If  the  signals  bandwidth  exceeds 
the  sampling  frequency  then  the  higher  frequency  components  overlap  and  sum  to 
introduce  errors  in  the  frequency  spectrum  of  the  time  signal.  This  is  known  as  aliasing. 
Figures  1  though  3  illustrate  the  concept  of  aliasing  and  the  Nyquist  sampling  rate. 

When  dealing  with  dynamic  real  time  signals,  the  bandwidth  of  the  signal  may  be 
quite  large.  In  order  to  effectively  sample  the  signal,  an  arbitrary  frequency  limit  must  be 
set  by  the  programmer.  Frequency  components  above  this  limit  are  usually  not  required  for 
systems  analysis  and  may  be  discarded  prior  to  sampling.  One  of  the  easiest  methods  to  do 
this  is  by  an  analog  low  pass  filter.  Setting  the  low  pass  filter  at  the  desired  frequency  limit 
allows  the  higher  frequency  components  of  the  signal  to  be  ignored.  Because  the  cut  on 
frequency  of  a  filter  is  defined  as  its  half  power  or  -3  db  point,  the  sampling  rate  is  usually 
somewhat  greater  than  twice  the  cutoff  frequency  of  the  filter.  This  allows  a  reduced 
chance  of  aliasing  the  signal  and  optimum  resolution  in  the  frequency  analysis. 

Another  method  of  avoiding  aliasing  is  to  sample  the  signal  at  three  to  four  times  its 
half  power  frequency  component.  This  ensures  that  any  aliasing  that  takes  place  is  in  the 
region  of  low  power  signals  and  that  it  does  not  affect  the  important  frequency  components 
of  the  signal.  The  major  drawback  to  this  method  is  that  by  sampling  at  a  higher  rate  the 
portion  time  signal  rate  is  decreased.  This  influences  resolution  of  the  spectrum. 


Digitized  Signet 

Figure  1 :  Sampled  Analog  Signal  in  the  Time  Domain 
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t  0  r 

Area  of  Aliasing 


Figure  3:  Signal  Sampled  at  less  than  the  Nyquist  Sampling  Frequency 


When  the  signal  is  transformed  into  its  frequency  components  via  a  discrete  FFT,  the 

minimum  frequency  resolution  is  inversely  proportional  to  the  time  record  length  of  the 
sampled  signal.  Thus  if  the  sampling  rate  is  Fs,  and  the  time  record  consist  of  n  samples, 

then  the  frequency  resolution  will  be: 


min 


1  Fs 

n  A  t  n 


[Eqn.  4] 


where  n  is  the  number  of  samples  in  the  data  set  The  n  A  t  term  is  the  time  record  length 
of  the  time  signal  and  is  noted  by  T. 
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Thus,  in  order  to  achieve  the  best  frequency  resolution,  the  longest  possible  time 
record  length  is  desired.  To  do  this  with  a  fixed  number  of  samples,  the  slowest  sampling 
frequency  possible  must  be  used.  A  filtered  signal  in  conjunction  with  a  slow  sampling 
rate  will  yield  non-aliased,  high  resolution  frequency  information. 

Fast  Fourier  transform  techniques  assume  that  the  time  signal  is  periodic  in  nature.  In 
order  to  meet  this  periodicity  requirement,  the  time  signal  must  be  weighted,  or  windowed. 
This  involves  tapering  the  information  by  various  methods  without  altering  the  frequency 
components  of  the  signal.  One  of  the  primary  windows  used  is  the  Hanning  or  cosine 
window.  This  multiplies  the  time  signal  by  a  sine  wave  of  amplitude  V2/3  and  a  period  of 
2T.  The  time  signal  is  then  zero  in  the  region  of  t  =  0  and  t  =  T.  This  introduces  a 
frequency  component  at  the  frequency  of  half  of  the  minimum  frequency  resolved  by  the 
Fourier  transform  method  which  will  not  affect  the  results  obtained  by  the  FFT.  Thus  the 
periodicity  requirements  are  met  and  the  frequency  components  of  the  time  signal  are  not 
altered. 

Important  additional  measurements  performed  on  the  mechanical  signatures  include 
the  power  spectra  density  function  (or  autospectra),  the  cross  spectra  density  function,  the 
system  transfer  function,  and  coherence.  The  input  power  spectra,  Gxx(co),  is  a  measure  of 
of  the  energy  content  of  the  input  signal  at  a  frequency  w.  Likewise,  the  output  power 
spectra,  Gyy(co),  is  a  measure  of  the  energy  content  of  the  output  signal.  Both  the  input 
power  spectra  and  output  power  spectra  are  real  valued  functions.  The  cross  power  spectra 
Gxy(co),  is  a  complex  value, 

GXy(<o)  =  GxyR(o>)  +  jGxyi(co),  [Eqn.  5] 


where  GxyR(«)  and  Gxyi(co)  are  the  real  and  imaginary  portions  of  the  cross  power 
spectra.  The  cross  power  spectra  is  a  measure  of  the  amount  of  energy  transferred  from  the 
input  point,  x,  to  the  output  point,  y,  at  a  frequency  to. 

The  system  transfer  function,  G(co),  is  the  ratio  of  the  Fourier  transform  of  the  input 


and  output  time  signals, 

r{y(0}  Y(o» 


[Eqn.  6] 


It  too  is  a  complex  number.  It  consist  of 
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G((o)  =  Gr  (w)  +  j  G  i(co), 


[Hqn.  7] 


where  Gr(co)  and  Gi((o)  are  the  real  and  imaginary  portions  of  the  transfer  function.  One 
method  of  calculating  the  transfer  function  uses  the  auto  spectra  and  cross  spectra 
information. 


Ig«o>|2  -  SodSSi 

Gxx(co) 


Gr(CD)  = 


GXVr(£0) 


[Eqn.  8] 


Gxx(ca) 

Gi(fl»  =  “^7“  •  [Eqn.  8] 

Gxx(co) 

The  coherence,  Yxy.  is  a  measure  of  the  noise  contamination  of  an  output  signal  y(t) 
from  sources  other  than  the  input  signal  x(t).  The  coherence  function  can  be  interpreted  as 
the  fractional  portion  of  the  output  spectrum  of  y(t)  which  is  linearly  due  to  x(t)  at 
frequency  to  [Ref.  3].  A  coherence  of  1  indicates  perfect  transmission  from  the  input  point, 
x,  to  the  output  point,  y.  This  indicates  that  there  is  no  noise  contamination.  A  coherence 
of  0  (zero)  indicates  severe  contamination  by  noise  at  the  output  point.  This  indicates  that 
none  of  the  signal  at  y  is  due  to  excitation  at  point  x.  The  coherence  is  calculated  by 


V^yfco)  = 


lGxy(C0)l2 
Gxx(CO)  Gyy(CO) 


[Eqn.  9] 


B.  THE  HILBERT  TRANSFORM  AND  THE  FREQUENCY  RESPONSE 

The  Hilbert  transform  has  been  identified  as  a  means  of  determining  the  degree  of 
linearity  of  a  system's  frequency  response.  It  has  been  shown  that  accurate  identification 
of  a  range  of  commonly  occurring  structural  non-linearities  such  as  friction,  stiffness  and 
power-law  damping  is  possible  in  systems  where  the  modes  are  well  separated  [Ref.  4]. 
By  taking  the  Hilbert  transform  of  a  systems  frequency  response  and  comparing  it  to  the 
original  frequency  response,  the  linearity  of  the  system  can  be  readily  determined.  For 
linear  systems  the  Hilbert  transform  of  the  frequency  response  H(co)  will  be  identical  to  the 
frequency  response  of  the  system  G(co). 

The  Hilbert  transform  of  a  real-valued  function  g(t)  is  defined  by: 


_i  fg( 

"*J  r- 


[Eqn.  10] 


In  the  frequency  domain  it  is  defined  by: 


H(co)  =  PV  f-^-dS, 

J  co-c 


[Eqn.  11] 


where  H(cd)  is  the  Hilbert  transform  of  the  complex  valued  G(£)  and  PV  is  defined  as  the 
Cauchy  Principle  Value  for  the  integral. 


G(£)  =  Gr(^)  +j  Gi(^) 


[Eqn.  12] 


where  Gr(^)  and  GI(^)  are  the  real  and  imaginary  parts  of  the  G<£).  It  is  important  to  note 
that  the  Hilbert  transform  is  in  the  same  domain  as  the  original  function.  Thus  the  Hilbert 
transform  of  a  real  valued  time  function  is  in  the  time  domain,  and  the  Hilbert  transform  of 

a  function  in  the  frequency  domain  is  in  the  frequency  domain.  The  Hilbert  transform 
H(oj)  is  complex -valued  in  the  frequency  domain. 


H(co)  =  Hr(g»  +  j  Hi(co), 

where  Hr(gj)  and  Hi(co)  are  the  real  and  imaginary  parts  of  H(co). 


[Eqn.  13] 


In  order  to  avoid  the  singularity  point  at  %  =  ©,  the  Hilbert  transform  the  integration 
must  be  divided  into  two  portions  with  limits  as  follows: 


—  Ur  C  oq 


co+  e 


Because  Gr(^)  is  an  even  function  and  Gi(^)  an  odd  function, 


GK)  =  G*($), 


[Eqn.  15] 
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where  *  denotes  the  complex  conjugate  of  the  complex  valued  function  G(£),  Equation  14 
must  be  further  subdivided  into  its  real  and  imaginary  parts: 


Hr(g>)  =  -  —  PV  f4y— t  d; 
*  </  5  “ 


[Eqn.  16] 


One  method  utilized  to  calculate  the  Hilbert  transform  involves  a  numerical 
integration.  Due  to  the  finite  frequency  resolution  of  a  discrete  Fourier  transform,  a 
numerical  integration  of  the  frequency  response  will  always  yield  inaccuracies  at  the 
singularity  point  to  =  £;.  The  equations  for  evaluating  the  Hilbert  transform  numerically  are 
given  by  [Ref.  5]: 


where  H(co)  denotes  the  Hilbert  transform  and  G(co)  is  the  system's  transfer  function.  The 
R  I 

correction  terms,  Ej  and  E  ■ ,  are  dependent  upon  the  system's  frequency  response  and  the 


range  of  the  frequency  response  under  examination. 

A  correction  term  is  needed  to  correct  for  resonant  peaks  outside  the  range  under 
examination,  and  others  to  correct  for  the  truncation  of  the  frequency  response  information 
at  the  limits  of  the  frequency  band  under  examination.  Corrections  for  multimode  systems 
are  usually  limited  to  a  shift  of  the  imaginary  portion  of  the  Hilbert  transform.  The 
estimate  of  the  shift  necessary  to  correct  this  is  given  by  [Ref.  6]: 
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[Eqn.  18] 


c,s  = 


j  w  X, 


where  CIs  is  the  correction  to  the  imaginary  portion  of  the  Hilbert  transform  due  to  the 
resonant  frequency,  tos.  Xs  is  the  the  imaginary  portion  of  the  modal  constant  for  the  s^ 
resonant  frequency.  An  estimate  of  the  value  of  X$  may  be  evaluated  from  the  the  fact  that: 

HiK)=^  [Eqn.  19] 

5stos 


Xs  =  H,(tos)  6S  cos. 


[Eqn.  20] 


8S,  the  structural  damping  ratio,  may  be  estimated  by  using  the  half  power  bandwidth  of  the 
resonant  frequencies  outside  the  band  under  investigation.  The  value  of  Hj(cos)  may  be 

taken  directly  from  the  imaginary  portion  of  the  transfer  function.  Thus  as  the  resonant 
frequencies  outside  the  band  under  investigation,  tos,  get  larger,  their  influence  of  it 
correction  term  CIs  diminishes.  For  structures  with  many  resonant  modes,  only  those 
modes  near  the  frequency  band  under  investigation  need  be  used  as  correction  terms. 

A  second  method  of  evaluating  the  Hilbert  transform  applied  in  this  study  is  the  use 
of  the  Fourier  transform  to  perform  the  Hilbert  transform.  The  integral 


LfS 

K  J  C 


[Eqn.  21] 


is  the  convolution  integral  in  the  frequency  domain.  Since  convolution  in  the  frequency 

domain  is  equivalent  to  multiplication  in  the  time  domain,  this  is  equivalent  to  taking  the 
Fourier  transform  of  G(co)  and  multiplying  it  by  the  Fourier  transform  of - . 

7t  tO 

The  Fourier  transform  of - is: 

71  0) 


1  1  f  *j  t>0 

-  I  =  -  j  sgn(t)  =  <  0  t=0  , 

7t  to  j  l  j  t<0 


The  Fourier  transform  of  G(^)  is: 


T{G($)j  =g(t). 


[Eqn.  22] 


[Eqn.  23] 
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There  exist  many  fast  Fourier  transform  routines  for  complex  vectors  in  computer 
libraries.  Thus  an  easier  method  to  calculate  the  Hilbert  transform  which  avoids  the 
numerical  integration  required  of  the  Cauchy  Method  is  to  take  the  Fourier  transform  of  a 
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2  t  >0 

real  signal  G(co)  and  multiply  it  by-!  1  t  =0  .  Then  by  taking  the  inverse  Fourier 

l  0  t  <0 

transform  of  this  modified  step  response,  the  transfer  function  G(co)  and  its  Hilbert 
transform  H(oj)  can  be  separated  out  [Ref  7].  Figure  4  summarizes  the  method  used  to 
calculate  the  Hilbert  transform  using  the  Fourier  transform,  Equations  21  through  27. 


Figure  4:  Fast  Hilbert  Transform  Flow  Diagram 

When  using  the  FFT  to  calculate  the  Hilbert  transform  a  FFT  routine  capable  of 
performing  inverse  Fourier  transforms  must  be  used.  Because  the  FFT  routine  assumes 

convolution  in  the  frequency  domain  has  taken  place,  the  Fourier  transform  is  symmetric 

T 

about  the  point  y,  where  T  denotes  the  time  record  length  of  the  input  signal.  The  inverse 

Fourier  transform  also  assumes  that  convolution  has  taken  place.  An  equivalent  method  of 

T 

setting  all  negative  frequency  component  to  zero  is  by  setting  all  components  above  y equal 

T 

to  zero.  All  component  from  f>0  toy  are  multiplied  by  2.  The  inverse  Fourier  transform 
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IS  then  performed  yielding  the  original  function  and  its  Hilbert  transform.  This  method  is 

faster  and  just  as  accurate  as  the  numerical  integration  method  for  systems  when  applied 

properly.  Correction  terms  are  still  required  when  using  the  FFT  method  of  calculating  the 
Hilbert  transform. 
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III.  EXPERIMENTAL  PROCEDURE 

A.  EQUIPMENT  UTILIZED 
1.  The  MASSCOMP  5400 

The  MASSCOMP  5400  computer  used  to  acquire  the  data  and  process  it  for  the 
experimental  work  was  configured  as  follows: 

1  71  Mbyte  hard  disc  drive 
1  300  kbyte  "floppy"  disc  drive 
6  Mbyte  Random  Access  Memory 

1  16  channel,  1  MHz  maximum  aggregate  sampling  rate  12  bit 
resolution  Analog  to  Digital  Converter  [AD12FA] 

1  16  channel,  333  kHz  maximum  aggregate  sampling  rate  12  bit 
resolution  Analog  to  Digital  Converter  [EF12M] 

5  Timing  Docks,  maximum  frequency  6  MHz 
1  (2)  channel  Digital  to  Analog  Converter  [EF12M] 

1  16  channel  sample  and  hold  temporary  storage  card  [SH16A] 

Various  display  and  output  devices 


The  capabilities  of  the  various  data  acquisition  devices  installed  in  the  computer  are 
summarized  in  Appendix  A. 

The  MASSCOMP  5400  is  a  Motorola  68020  based  32  bit  word  length  system. 
It  is  run  on  the  UNIX  operating  system.  Language  capabilities  include  "C"  and 
FORTRAN.  Software  use  to  operate  the  system  includes  a  graphics  presentation  package, 
Dynamic  Signal  Analysis  Subroutines,  Window  Presentation  Packages  and  the  software 
necessary  to  operate  the  Data  Acquisition  devices. 

The  data  analysis  routines,  except  where  noted,  were  copied  from  those 
provided  by  the  MASSCOMP  Dynamic  Signal  Analysis  Package  (SP-60).  These  were 
originally  copied  from  the  collection  of  Programs  for  Digital  Signal  Processing  from  the 
IEEE  Acoustics,  Speech  and  Signal  Processing  Society. 

The  programming  necessary  to  collect  and  process  the  data  was  written  in 
FORTRAN.  All  programs  referred  to  in  the  remainder  of  this  paper  are  documented  in  the 
appendices.  It  should  be  noted  that  the  software  packages  for  the  data  acquisition  and 
presentation  portions  include  subroutine  calls  that  are  unique  to  the  MASSCOMP  system 
and  are  written  specifically  to  support  MASSCOMP  hardware.  Further  investigation  into 
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system 
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the  software  documentation  provided  by  MASSCOMP  is  warranted  for  a  complete 
understanding  of  the  data  collection  and  presentation  routines. 

In  addition  to  the  MASSCOMP  specific  equipment,  various  recording  devices, 
dynamic  signal  analyzers  and  signal  generators  were  used  to  provide  the  necessary  data  and 
a  means  of  verification  of  the  accuracy  and  applicability  of  the  MASSCOMP  5400  system. 

2.  Analog  to  Digital  Conveners 

Analog  to  Digital  converters  work  in  several  different  manners.  While  the 

method  used  to  digitize  a  signal  is  not  important,  the  ultimate  goal  of  a  converter  is  to 

compare  a  sampled  voltage  and  to  digitize  it  to  the  nearest  voltage  value  within  its  dynamic 

range.  A  twelve  bit  A/D  converter  is  able  to  separate  its  dynamic  range  into  212  individual 
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levels.  Thus  a  converter  with  a  20  volt  dynamic  range  is  able  to  obtain  a  resolution  of 

or  0.00488  volts.  A  12  bit  converter  with  only  a  10  volt  range  will  have  twice  this 
resolution.  The  speed  of  the  conversion  is  dependent  upon  the  method  of  digitization;  the 
fastest  devices  having  speed  on  the  order  of  20  nanoseconds  and  the  slower  ones  of  1.2 
microseconds  [Ref.  8  p.  97]. 

The  transfer  rate  of  the  data  to  a  usable  state  is  dependent  upon  the 
programming  of  the  converter.  There  are  three  methods  available  to  the  programmer  for 
data  acquisition.  The  first  is  reading  the  information  directly  into  an  array  for  later  use  in 
the  program.  This  is  the  easiest  and  usually  the  quickest.  If  the  data  is  read  into  dynamic 
memory  then  the  converter  can  be  run  at  its  maximum  speed  or  at  the  data  bus’  speed 
whichever  is  less.  Hie  second  method  is  a  direct  disc  transfer.  This  is  usually  slower  than 
the  array  transfer  method.  If  the  information  is  to  be  read  into  a  hard  disc  for  later  retrieval 
and  reduction  then  the  programmer  is  usually  limited  to  the  transfer  rate  capability  of  the 
particular  hard  disc  system  available.  One  drawback  of  the  direct  to  disc  transfer  is  that  the 
programmer  has  no  access  to  the  information  until  after  the  transfer  is  complete.  The  third 
and  more  complicated  method  is  to  read  the  data  into  memory  buffets.  The  advantage  of 
this  method  is  that  the  data  can  be  manipulated  as  the  buffers  are  filled  rather  than  having  to 
wait  until  the  final  array  or  disc  transfer  is  complete.  Most  micro  computers  allow  a 
programmer  to  select  any  one  of  the  three  possibilities  depending  upon  the  application  of 
the  data. 

Once  the  data  was  transferred  to  the  hard  disc  in  either  the  raw  or  reduced  form 
it  could  be  transferred  to  magnetic  tape  or  floppy  disc.  This  allows  the  volume  of  the 
information  on  the  hard  disc  to  be  substantially  reduced  after  the  acquisition  and  analysis  is 
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completed.  If  the  data  is  required  at  a  later  date  it  is  available  for  transfer  from  the 
permanent  recordings  (floppy  disc  or  tape). 

3.  Digital  to  Analog  Conveners 

Digital  to  Analog  converters  work  in  somewhat  the  same  manner  of  the  A/D 
converter,  only  their  task  is  simpler.  Because  the  goal  is  to  send  out  an  analog  voltage, 
there  is  no  need  to  compare  the  desired  value  to  the  output  signal.  Because  the  D/A 
convener  is  essentially  a  resistor  capacitor  circuit,  the  settling  time  requirement  of  the  RC 
circuit,  and  thus  the  D/A  converter  determines  its  speed  capability.  Settling  times  are 
dependent  upon  the  exact  circuitry  and  the  resolution  of  the  converter.  Eight  bit  converters 
may  have  settling  times  of  0. 1  microseconds,  while  sixteen  bit  converters  may  have  settling 
times  of  75  microseconds  [Ref.  8  p.  95],  The  dynamic  range  of  D/A  converters  are 
comparable  to  that  of  A/D  converters.  The  programmer  may  also  use  the  same  three 
methods  of  transfer  (array,  disc  and  buffer)  for  digital  to  analog  conversions. 

Analog  to  Digital  converters  and  Digital  to  Analog  converters  require  timing 
clock  signals  in  order  to  determine  the  frequency  of  sampling  or  output.  Programmable 
clock  chips  are  used  to  provide  voltage  pulse  trains  for  the  devices.  Thus  the  programer  is 
able  to  set  the  sampling  rate  of  an  A/D  device  and  the  output  speed  of  a  D/A  device.  With 
clock  chips  capable  of  providing  6  MHz  pulse  trains,  the  programmer  is  limited  only  by  the 
time  characteristics  of  his  conversion  devices  and  the  computers  systems  data  bus  capacity. 

4.  Verification  Procedure 

The  process  of  determining  whether  the  MASSCOMP  5400  system  was 
suitable  for  the  near  real  time  processing  of  shock  and  vibration  signals  was  completed  in  a 
series  of  steps.  First,  using  analog  data  pre-recorded  on  magnetic  tape  from  underwater 
explosions,  the  analog  converters  were  tested  for  accuracy  and  speed.  The  MASSCOMP 
sampled,  recorded,  analyzed  and  displayed  the  information.  Comparisons  to  the  data 
analyzed  by  other  means  were  used  to  verify  the  accuracy  of  the  acquired  data.  Once 
validation  of  sampling  speed  and  analysis  of  selected  data  sets  was  accomplished  all  data 
sets  were  transferred  to  the  hard  disc  for  archival  purposes. 

The  second  method  of  systems  verification  was  completed  by  impact  testing  of 
an  aluminum  plate.  Data  analysis  was  expanded  to  include  power  spectral  analysis  and 
coherence  measurements.  Again  the  results  were  compared  to  those  obtained  by  an 
separate  HP  3562A  Dynamic  Signal  Analyzer  in  order  to  verify  the  accuracy  of  the 
MASSCOMP  system. 
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5.  Application  of  the  MASSCOMP  5400  System 

Once  the  system  was  determined  to  operate  satisfactorily,  it  was  used  on  the 
vibrational  analysis  of  a  linear  response  system.  In  addition  to  the  frequency  response, 
power  spectral  density  and  coherence  measurements,  the  Hilbert  transform  was  used  to 
examine  the  degree  of  linearity  of  the  frequency  response.  This  enabled  a  rapid 
identification  of  any  nonlinearities  in  the  system  transfer  function.  The  two  methods  of 
computing  the  Hilbert  transform,  numerical  integration  and  use  of  the  Fourier  transform, 
were  compared  for  accuracy.  Once  the  Hilbert  transform  method  was  chosen,  the 
frequency  response  characteristics  of  a  non-linear  system  were  evaluated  to  ensure  the 
Hilbert  transform  identified  the  non-linearities. 

B.  SHOCK  AND  VIBRATION  ANALYSIS 

1  ■  Analysis  of  Prerecorded  Underwater  Explosion  Data 

The  analysis  of  prerecorded  underwater  explosion  data  involved  the  transfer  of 
pressure  and  strain  gage  records  from  a  Honeywell  Model  101  sixteen  channel  FM  tape 
recorder  to  the  MASSCOMP  system.  Data  was  recorded  at  a  tape  speed  of  120  inches  per 
second.  The  AD12FA  Analog  To  Digital  convener  card  was  used  for  the  transfer  of  data. 
Figure  5  shows  the  schematic  diagram  of  the  connections  used  for  the  data  transfer. 

The  bandwidth  capability  of  the  tape  recorder  was  250  kHz.  The  bandwidth  of 
the  pressure  transducer  data  was  250  kHz,  while  the  strain  gage  data  was  only  10  kHz. 
Location  of  data  on  the  tape  was  noted  by  a  trigger  signal  on  one  of  the  sixteen  data 
channels.  During  the  explosive  testing  procedure,  an  electrical  wire  with  a  DC  voltage 
applied  to  it  was  wrapped  around  the  explosive  charge.  As  the  charge  was  detonated,  the 
cable  was  severed  and  the  voltage  sensed  at  the  recorder  dropped  to  zero.  Monitoring  the 
trigger  channel  for  this  drop  in  the  voltage  enabled  the  program  to  start  the  clock  module. 
When  the  clock  began  sending  pulses  to  the  Analog  To  Digital  converter,  information  was 
transferred  off  the  tape  into  the  data  storage  arrays  in  the  computer  program. 

The  Analog  to  Digital  Convener  section  of  the  EF12M  multifunction  card  is 
limited  to  a  bandwidth  of  only  160  kHz.  This  required  the  use  of  the  AD12FA  Analog  To 
Digital  convener.  This  allowed  sampling  of  the  recorded  data  at  a  maximum  frequency  of 
1.0  MHz,  well  beyond  the  bandwidth  of  the  original  signal.  Figure  6  is  the  pressure 
transducer  data  set  used  as  a  reference  to  compare  the  data  recorded  by  the  A/D  system. 
This  was  recorded  by  a  HP5451C  computer  system  at  a  sampling  rate  of  500  kHz.  Data 
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was  acquired  on  the  MASSCOMP  system  at  various  sampling  rates  as  shown  in  Table  1 
and  then  compared  to  the  reference  record.  The  program  utilized  to  acquire  the  data  is 
shown  in  Appendix  B. 

TABLE  1: 

SUMMARY  OF  SAMPLING  RATES  USED  TO  EXAMINE  UNDERWATER 

EXPLOSION  DATA 


Tape  Speed  (ips) 

Sampling  Rate 

120 

1  MHz 

120 

750  kHz 

120 

500  kHz 

120 

250  kHz 

Honeywell  Model 
101  FM 

Tape  Recorder 
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Figure  6:  Reference  Pressure  Record  Used  as  the  Standard 

Figures  7  through  10  are  the  data  recordings  from  the  same  pressure  pulse 
recorded  at  the  sampling  rates  listed  above.  Comparisons  of  these  Figures  to  the  base  line 
record  revealed  that  the  500  kHz  sampling  rate  matched  the  reference  ideally.  This  was  as 
predicted  by  the  Nyquist  sampling  theorem.  Recordings  at  rates  less  than  the  Nyquist 
frequency  missed  portions  of  the  data  that  was  available  on  the  tape.  Rates  faster  than  the 
Nyquist  frequency  exceeded  the  bandwidth  of  the  tape  recorder  and  the  reference  record 
and  subsequently  appeared  "noisy".  Because  there  was  no  method  to  determine  if  this 
noise  was  actually  recorded  information  or  noise  produced  by  the  tape  recording  and  play 
back  process,  the  sampling  rate  for  all  further  analysis  of  the  prerecorded  underwater 
explosion  data  was  fixed  at  500  kHz.  The  strain  gage  data  was  then  acquired  from  the 
Honeywell  tape  recorder  and  transferred  to  the  MASSCOMP  hard  disc. 

The  data  was  transferred  to  the  the  MASSCOMP  system  one  track  at  a  time.  In 
order  to  utilize  the  system,  as  it  is  presently  configured,  for  real  time  data  acquisition  in  the 
field  the  sixteen  channel  sample  and  hold  card  (SH16F)  and  the  single  function  Analog  to 
Digital  Convener  (AD12FA)  must  be  used.  The  only  limitation  is  the  data  transfer  bus, 
which  is  limited  to  2  Megabytes  per  second  This  allows  one  million  sixteen  bit  words  per 
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second  to  be  gathered  in  real  time.  Table  2  summarizes  the  bandwidth  limitations  based 
upon  the  number  of  sampled  channels. 

TABLE  2:  BANDWIDTH  LIMITATIONS  FOR  MULTIPLE  DATA  CHANNELS 


Number  of  Channels  Sampled 

Bandwidth  Limit 

1 

500  kHz 

2 

250  kHz 

4 

125  kHz 

8 

62.5  kHz 

16 

31.5  kHz 

Thus  it  would  be  impossible  to  record  more  than  two  pressure  pulses  in  real 
time  at  a  500  kHz  sampling  rate.  It  would  be  possible  to  record  up  to  sixteen  strain  gage 
traces  at  a  20  kHz  sampling  rate  on  the  AD12FA  A/D.  The  present  system  configuration 
limits  the  data  acquisition  to  one  type  of  data,  either  pressure  or  strain  gage.  The  single 
sample  and  hold  card  can  only  be  clocked  at  one  frequency,  limiting  the  choice  to  one  type 
of  data.  A  second  sample  and  hold  card,  connected  to  the  EF12M  multifunction  module, 
would  allow  the  simultaneous  sampling  of  both  pressure  and  strain  gage  information. 

In  addition,  to  utilize  the  system  for  live  data  recording  the  program  listed  in 
Appendix  B  would  have  to  be  modified  to  include  the  use  of  the  sample  and  hold  card. 
This  would  necessitate  a  modification  to  the  identification  of  the  channels  to  be  sampled  and 
the  clocking  method  used.  Due  to  the  increase  in  the  volume  of  data,  the  storage  arrays 
would  have  to  be  modified  to  allow  up  to  sixteen  channels  of  data.  This  would  create  the 
need  for  large  arrays  (16  [one  for  each  channel]  by  20,000  [one  second  of  data  at  20  kHz 
sampling  rate])  for  the  raw  data.  Another  method  would  be  to  read  the  information  directly 
to  the  disc  and  then  process  it  after  the  data  is  on  the  disc.  This  would  still  allow 
processing  of  the  data  in  the  field  and  allow  analysis  within  minutes  of  the  explosion. 


2.  Impact  Testing  of  an  Aluminum  Plate 

Using  the  equipment  schematic  diagram  shown  in  Figure  1 1,  the  digital  signal 
analysis  portion  of  the  data  acquisition  system  was  tested.  A  modally  tuned  impact  hammer 
with  a  bandwidth  of  1,200  Hz  was  used  to  strike  the  aluminium  plate.  A  force  transducer 
located  in  the  hammer  provided  a  voltage  signal  to  one  channel  of  the  A/D  converter.  An 
accelerometer  located  on  the  plate  provided  the  plate  acceleration  data  to  a  second  channel  of 
the  A/D  converter. 

The  data  was  sampled  at  a  rate  of  5,000  Hz  per  channel  providing  a  2,500  Hz 
bandwidth  for  the  sampled  signal.  Because  the  energy  level  of  the  input  signal  was 
negligible  beyond  2,500  Hz,  a  low  pass  filter  was  not  required.  Data  blocks  of  1024 
samples  were  used  yielding  a  frequency  resolution  of  4.88  Hz.  This  ensured  that  an 
adequate  number  of  data  points  were  in  the  frequency  band  width  of  the  hammer.  The 
program  utilized  to  acquire  the  data  is  shown  in  Appendix  C. 

The  time  signals  for  the  hammer  and  plate  were  recorded  and  averaged  in  the 
frequency  domain.  The  frequency  information  was  derived  utilizing  the  FFT  routine 
FFT842  [Ref.  9],  while  the  power  spectral  information  and  coherence  measurements  were 
derived  by  utilizing  the  CCSE  routine  [Ref.  10].  Both  of  these  routines  were  modified  to 
suit  the  particular  characteristics  of  the  impact  testing  procedure  and  the  microcomputing 
system. 

The  dynamic  signal  analysis  was  first  performed  on  a  single  data  record.  The 
coherence,  transfer  function  and  frequency  components  of  the  signals  were  computed.  As 
expected  the  coherence  for  a  single  impact  was  1  throughout  the  frequency  band  under 
examination.  The  identical  test  was  recorded  using  a  HP3562A  Dynamic  Signal  Analyzer 
for  comparison.  The  analysis  performed  utilizing  the  super  microcomputer  yielded 
comparable  results.  Figures  12  through  14  display  the  recorded  data  for  this  single  impact 
test  from  both  sources. 
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Figure  1 1 :  Experimental  Set  Up  for  Impact  Testing  of  Aluminum  Plate 


The  impact  testing  was  then  averaged  to  determine  the  influence  of  averaging  on 
the  transfer  function.  Figures  15  through  18  are  the  data  recordings  for  sixteen  impacts. 
As  expected  the  frequency  components  of  both  the  hammer  and  accelerometer  data 
smoothed  considerably  when  averaged.  The  coherence  decreased  in  areas  of  noise 
contamination  and  non-linearities. 


Steady  state  vibrations  were  measured  using  the  apparatus  shown  in  the 
schematic  diagram  in  Figure  19.  An  aluminum  beam  was  mounted  so  that  it  could  be 
excited  by  a  shaker.  The  input  acceleration  was  measured  at  the  mounting  jaws  and  the 
output  acceleration  was  measured  at  the  tip  of  the  beam.  The  beam  was  excited  using  a 
random  noise  excitation  of  1  kHz  and  500  Hz  bandwidths.  The  system  transfer  function 
was  calculated,  and  the  Hilbert  transform  of  the  real  and  imaginary  portions  of  the 


frequency  response  was  performed  to  determine  if  the  system  behaved  in  a  linear  manner. 
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Figure  12:  Input  Power  Spectra  for  Single  Impact 
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Figure  13:  Output  Power  Spectra  for  Single  Impact 
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Figure  18:  Coherence  Measurement  for  Multiple  Impact 
A  low  pass  Filter  ensured  that  no  frequency  components  above  the  range  of 
investigation  contaminated  the  transfer  function.  The  sampling  rate  was  such  that  a 
frequency  resolution  of  1.3  Hz  was  obtained.  For  steady  state  vibration  response 
measurements,  a  modified  version  of  the  program  used  for  impact  testing  was  used.  The 
program  is  listed  in  Appendix  D.  Because  the  same  data  reduction  and  analysis  was 
performed  on  the  steady  state  vibrations  testing  as  was  done  on  the  impact  testing 
information  the  same  reduction  routines  were  used 

Figure  20  is  the  input  power  spectra  for  a  random  noise  excitation  at  a  1  kHz 
bandwidth  centered  at  500  Hz.  Figures  21  and  22  are  the  output  power  spectra  and 
coherence  measurements  for  the  same  excitation.  The  transfer  function  for  this  excitation  is 
shown  in  Figure  23.  The  use  of  random  noise  produced  a  low  coherence  measurement  in 
the  region  above  200  Hz.  In  order  to  improve  the  coherence  measurements  and  smooth  the 
transfer  function  in  the  low  frequency  area,  a  random  noise  excitation  of  only  500  Hz 
bandwidth  starting  at  zero  was  used.  The  results  of  using  this  excitation  are  shown  in 
Figures  24  through  29. 

An  attempt  to  utilize  a  swept  sine  wave  as  the  excitation  source  was  made.  Two 
methods,  one  using  a  HP  3562A  Dynamic  Signal  Analyzer  and  a  second  using  a  voltage 
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controlled  oscillator  driven  by  the  digital  to  analog  convener  of  the  EF12M  multifunction 
module,  were  used  to  sweep  through  a  frequency  band.  Using  the  Hewlett  Packard 
dynamic  signal  analyzer  was  not  effective.  As  the  dynamic  signal  analyzer  swept  through 
the  frequency  bandwidth,  the  MASSCOMP  would  sample  the  responses.  However,  due  to 
the  length  of  time  required  by  the  MASSCOMP  to  process  the  data,  the  MASSCOMP 
would  miss  some  frequencies  as  the  HP  Dynamic  Signal  Analyzer  swept  through  them. 
The  method  of  using  a  voltage  controlled  oscillator  was  more  effective;  however,  the 
processing  time  was  such  that  only  a  200  Hz  bandwidth  could  be  swept  without  skipping 
some  frequencies.  This  required  approximately  twenty-five  minutes  to  sweep  this 
bandwidth,  while  with  a  Dynamic  Signal  Analyzer  the  entire  bandwidth  from  0  to  1  kHz 
could  be  sampled  in  the  same  amount  of  time.  This  necessitated  the  use  of  random  noise  as 
the  only  excitation  source  of  the  shaker. 

The  results  of  the  steady  state  analysis  were  compared  to  the  HP  3562A 
Dynamic  signal  Analyzer  for  accuracy.  Figures  30  through  34  are  the  results  of  those  test. 
When  compared  to  the  results  gained  via  the  microcomputer,  the  differences  are  quite 
noticeable.  Of  particular  importance  is  the  difference  in  magnitude  in  the  input  power 
spectra,  the  output  power  spectra  and  the  cross  spectra.  While  the  shape  of  the 
measurement  results  are  nearly  identical,  the  magnitudes  are  different.  This  leads  to  an 
incorrect  approximation  of  the  real  and  imaginary  portions  of  the  transfer  functions.  This  is 
evident  in  Figures  35  through  38.  Figures  35  and  36  are  the  real  and  imaginary  portions  of 
the  transfer  function  measured  with  the  microcomputing  system.  The  identical 
measurements  using  the  dynamic  signal  analyzer  are  shown  in  Figures  37  and  38.  The 
difference  in  magnitude  is  immediately  apparent. 


Low  Pass  Filter 


Input  Accelerometer 
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Figure  19:  Steady  State  Vibration  Measurement  Schematic  Diagram  (upper)  and 

Cantilevered 

Beam  used  in  Experimental  Procedure  (lower). 
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Figure  23:  Transfer  Function  for  1  kHz  Random  Noise  Excitation 
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Figure  26:  Coherence 
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Figure  28:  Real  Portion  of  Cross  Power  Spectra  for  500  Hz  Random  Noise  Excitation 
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Figure  29:  Imaginary  Portion  of  Cross  Power  Spectra  for  500  Hz  Random  Noise 
3  '  Excitation. 
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Figure  30:  Input  Power  Spectra  for  500  Hz  Excitation  Measured  with  HP  3562A 


Figure  31:  Output  Power  Spectra  for  500  Hz  Random  Excitation  Measured  with  HP 

3562A 
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Figure  33:  Real  Portion  of  Cross  Power  Spectra  for  500  Hz  Random  Excitation  Measured  with  HP  3562A 
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Figure  34:  Imaginary  Portion  of  Cross  Power  Spectra  for  500  Hz  Random  Excitation  Measured  with 

HP3562A 
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r:  Real  Portion  of  Transfer  Function  for  500  Hz  Excitation  Measured 


3562A 


The  aluminum  beam  used  for  steady  state  vibrations  analysis  was  also  used  to 
evaluate  the  feasibility  of  including  the  Hilbert  transform  as  a  method  of  identifying  system 
irregularities.  There  are  two  resonant  modes,  one  at  137  Hz  and  the  second  at  891  Hz, 
shown  previously  in  Figure  23.  Using  the  half  power  bandwidth  method  to  estimate  the 
viscous  damping  ratio  at  the  second  resonant  frequency: 


_1  f0>2  ~ 


[Eqn.  28] 


where  co2  and  o)j  are  the  half  power  frequencies  for  that  resonant  point  and  cos  is  the 

resonant  frequency.  Using  the  information  for  the  second  resonant  frequency  from  Figure 
39,  the  viscous  damping  ratio  £s  was  determined  to  be  0.00180.  This  damping  ratio  was 

then  used  to  estimate  the  correction  terms  necessary  for  the  Hilbert  transform.  Recalling 
that  the  structural  damping  ratio  is  proportional  to  the  viscous  damping  ratio, 


5s=y 


[Eqn.  29] 


and  substituting  into  Equation  20, 


Xs  =  H,(cos)5scos  =  HjfcOj)  y»s 


[Eqn.  30] 


Substituting  Equation  30  into  Equation  18,  the  correction  term  for  the  imaginary  portion  of 
the  Hilbert  transform  is  then: 


r  _ ■  ■  CO  H,(0)s)  (0)2—0),) 

Ci,(0))  J  2  2  J  2  2  ~  J  2  2\  [Eqn.  31] 

C0S  -  CO2  2  (co2  -  CO2)  4  (cos  -  CO2) 

Using  the  damping  information  from  the  second  resonant  frequency,  the  correction  terms 
were  less  than  10' 5  and  were  considered  negligible  and  were  not  included. 
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The  real  and  imaginary  portions  of  the  original  transfer  function  and  the  Hilbert 
transform  of  the  transfer  function  are  shown  in  Figures  40  through  44.  The  numerical 
integration  method  of  performing  the  Hilbert  transform  was  not  accurate  in  the  area  of  the 
first  resonant  frequency.  Figures  45  and  46  show  that  the  Hilbert  transform  of  the  transfer 
function  was  not  significantly  different  than  the  original  transfer  function.  The  numerical 
integration  method  will  always  yield  incorrect  answers  for  systems  of  low  damping  at 
resonance. 

The  FFT  method  of  performing  the  Hilbert  transform  was  also  inaccurate.  Due  to  the 
transform  into  the  time  domain  and  subsequent  truncation,  there  appeared  to  be  high 
frequency  components  introduced  to  the  transfer  function.  This  is  evident  in  Figures  47 
and  48.  However,  the  shape  of  both  the  real  and  imaginary  portions  of  the  Hilbert 
transform  are  identical  to  the  original  transfer  function. 

When  using  the  frequency  response  of  a  model  of  a  linear  system  the  Hilbert 
transform  was  able  to  indicate  that  the  system  behaved  in  a  linear  fashion.  However,  due 
to  the  inaccuracy  in  measuring  the  cross  power  spectra  the  transfer  functions  used  in 
performing  the  Hilbert  transform  are  inaccurate.  Thus  in  order  to  ensure  that  the  Hilbert 
transform  in  yielding  correct  information,  the  original  transfer  function  must  be  correct; 
otherwise,  the  Hilbert  transform  will  indicate  a  system  that  is  in  fact  a  linear  system  as  a 
non-linear  system. 
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Figure  4 1 :  Imaginary  Portion  of  the  Transfer  Function  of  the  First  Resonant  Frequency 
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Figure  42:  Hilbert  Transform  of  the  Real  Portion  of  the  Transfer  Function  using  Numerical  Integration 
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A.  UNDERWATER  EXPLOSION  TESTING 

1 .  The  super  microcomputer  is  able  to  record  up  to  two  channels  of  shock  pressure  with 
250  kHz  bandwidth  or  sixteen  channels  of  strain  information  with  10  kHz  bandwidth 
directly  to  data  arrays.  It  is  possible  to  reduce  and  display  the  data  prior  to  storage  to 
a  recording  media. 

2.  In  the  present  configuration  of  only  one  sample  and  hold  device  for  the  data 
acquisition  hardware  only  one  type  of  data,  shock  pressure  or  strain,  can  be  recorded. 
The  addition  of  a  second  sample  and  hold  card  would  allow  two  different  data  types 
to  be  recorded  and  analyzed. 

B .  IMPACT  TESTING  OF  AN  ALUMINUM  PLATE 

1 .  The  data  acquisition  capabilities,  tied  to  the  Digital  Signal  Processing  software,  allow 
for  a  complete  dynamic  signal  analysis  of  a  structure  using  a  microcomputer. 

2.  Sampling  rates  of  32  kHz  are  sufficient  for  multi-channel  low  frequency  vibrations 
analysis.  However,  the  requirement  to  use  an  analog  filter  prior  to  sampling  requires 
the  use  of  a  filter  for  each  data  channel  monitored.  This  would  require  a  rack  of 
filters  and  reduce  the  ease  of  movement  of  the  system.  Setting  up  of  a  vibrations 
analysis  station  on  a  permanent  basis  would  be  required. 

C.  STEADY  STATE  VIBRATION  ANALYSIS 

1 .  The  data  acquisition  techniques  are  sound;  however,  the  microcomputer  data 
acquisition  equipment  yields  information  that  is  significantly  different  in  magnitude 
when  compared  to  a  proven  measurement  device.  This  is  due  to  the  12  bit  resolution 
over  the  10  volt  dynamic  range  of  the  analog  to  digital  converter. 

2.  The  coherence  measurements  using  an  external  random  noise  source  indicate  that  this 
is  not  the  optimum  method  of  performing  response  measurements.  Attempts  to 
improve  the  coherence  measurements  by  using  a  voltage  controlled  oscillator 
controlled  by  the  computer  as  a  sinusoidal  excitation  source  were  unsuccessful. 

3.  The  identification  of  non-linear  response  using  the  Hilbert  transform  is  limited  to 
systems  with  a  few  widely  spaced  resonant  modes.  It  requires  data  reduction  to  be 
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performed  by  the  user  prior  to  attempting  the  identification  of  non-linearhies. 
Because  of  the  poor  voltage  resolution  of  the  analog  to  digital  converters,  the  system 
transfer  information  used  as  an  input  for  the  Hilbert  transform  is  in  error.  Thus  the 
reliability  of  the  Hilbert  transform  for  use  determination  of  the  linearity  of  the 
frequency  response  is  poor. 

4.  The  present  data  presentation  in  graphical  format  capability  is  rudimentary  and 
requires  improvement.  The  capability  to  include  cursor  input  into  the  program  for 
axes  scaling  and  data  zooming  and  recording  would  greatly  improve  the  analysis 
capabilities  available  to  the  user.  The  capabilities  of  using  a  "mouse''  or  a  digital 
graphics  pad  as  an  control  device  to  delineate  points  of  interest  for  further  reduction 
would  greatly  enhance  the  capability  to  process  information. 


A.  UNDERWATER  EXPLOSION  TESTING 

1 .  Include  the  use  of  a  second  sample  and  hold  card.  This  would  allow  simultaneous 
recording  of  a  single  high  frequency  shock  pressure  signal  and  sixteen  low  frequency 
strain  gage  signals. 

B  IMPACT  TESTLNG 

1 .  Improve  the  presentation  software  so  that  cursor  input  via  a  "mouse"  can  be  used  to 
select  areas  of  interest  for  further  analysis. 

C.  STEADY  STATE  VIBRATION  ANALYSIS 

1  Explore  the  possibility  of  using  a  proved  data  acquisition  device,  the  HP  3562A,  in 
conjunction  with  the  programming  capability  of  the  MASSCOMP  system.  The 
EF12M  Multifunction  Module  has  an  IEEE  488  interface  capability.  The  HP  3562A 
may  be  able  to  acquire  the  data  and  transfer  it  to  the  MASSCOMP  for  exploration  of 
the  degree  of  linearity  of  the  transfer  function  by  use  of  the  Hilbert  transform. 

2  Reduce  the  processing  time  requirements  for  large  data  arrays.  FORTRAN 
programming  is  not  the  best  method  to  process  information  in  real  time  applications. 
Transferring  the  code  to  assembly  language  or  machine  code  would  greatly  reduce  the 
time  requirements  of  data  analysis  and  allow  faster  structure  analysis. 

*  Improve  the  presentation  software  as  noted  above  for  impact  testing. 

4  Develop  a  method  of  swept  sine  wave  excitation  of  the  structure  using  the 
microcomputer  as  the  controller.  The  digital  to  analog  converter  allows  for  two 
channels  of  control  information  to  be  sent  from  the  computer,  yet  it  must  operate 
independent  of  the  analog  To  digital  converter. 

3  Add  a  signal  conditioner  to  boost  the  level  of  the  input  signals  to  the  analog  to  digital 
converter.  This  would  result  in  more  accurate  calculation  of  the  system  transfer 
function  and  increase  the  reliability  of  the  results  of  the  Hilbert  transform 


measurements. 


APPENDIX  A 


MASSCOMP  5400  HARDWARE  CAPABILITIES 


The  capabilities  of  the  data  acquisition  devices  installed  on  the  MASSCOMP  5400  and 
the  microcomputer  itself  are  summarized  below; 

MASSCOMP  5400  Super  Microcomputing  system 
Motorola  68020  CPU  module 
6  Mbyte  of  system  memory 
Data  Acquisition  Hardware 

AD12FA  Analog  to  Digital  Converter 
EF12M  Multifunction  Module 
SH16F  Sample  and  Hold  Card 
Two  Display  Terminals 
HP  plotter 
Dot  matrix  printer 

AD12FA  Analog  to  Digital  Converter 
[Path  identified  as  /dev/dacpO/adfl] 

16  channel  input  single  ended,  8  differential  [channels  0-15] 

Programmable-gain  amplifier  (with  gains  selectable  by  software  of  1,2,4,  or  8) 

Sample  and  Hold  circuit  for  connection  to  the  SH-16F 

Maximum  aggregate  sampling  rate:  1  Mhz 

Resolution:  12  bits  (zero  padding  for  total  sample  length  of  16  bits) 

Input  Voltages: 

0  to  +10  V 
-5  to  +  5  V 

EF12M  Multifunction  Module 

Analog  to  Digital  Converter  Section 
[Path  identified  as  /dev/dacpO/adfO] 

16  channel  input  single  ended,  8  differential 

Programmable-gain  amplifier  (with  gains  selectable  by  software  of  1,2,4,  or  8) 

Maximum  aggregate  sampling  rate:  333  kHz 

Resolution:  12  bits  (zero  padding  for  total  sample  length  of  16  bits) 

Input  Voltages: 

-10  to  +10V 


Digital  to  Analog  Converter  Section 
[  Path  identified  as  /dev/dacpO/dafO] 

2  channels 

Maximum  clock  rate  per  channel:  333  kHz 

Resolution:  12  bits  (zero  padding  for  total  sample  length  of  16  bits) 

Output  Voltages:  -10  to  +10  V 
Parallel  I/O  Section 

16  Channels  Input  and  Output  utilizing  hardware  handshake  for  transmission 
of  data 

Clock  Section 

f  Path  identified  as  /dev/dacpO/efclkn  (n  identifies  clock  0-5)] 

5  Counters  with  output,  source  input  and  gate  input  connections 
Frequency  capability  ~0  Hz  to  3  MHz 

SH-16F  Sample  and  Hold  Module 

[Path  /dev/dacpO/adfl,  channels  16  through  31  on  the  AD12FA  card] 

16  Single-ended  Input  Channels  or, 

8  Differential  Channels 

Maximum  of  3  SH-16F  modules  may  be  connected  to  a  AD12FA  or  EF12M  to 
give  the  capability  of  sampling  a  maximum  of  48  channels  with  a  single  A/D 
converter.  The  maximum  aggregate  sampling  rate  must  never  exceed  the 
capability  of  the  "host"  A/D  converter. 
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UNDERWATER  EXPLOSION  DATA  TRANSFER  PROGRAM 


This  appendix  contains  the  program  used  to  collect  the  data  of  the  underwater  explosion 
recordings  from  the  FM  tapes.  All  routines  were  written  by  the  author.  For  further 
explanation  of  the  data  acquisition  specific  routines  the  reader  is  encouraged  to  refer  to  the 
MASSCOMP  documentation  "  Data  Acquisition  User's  Manual."  Figure  46  is  the  program 
flow  diagram.  Sections  of  the  program  follow  the  general  outline  presented  in  the  flow 
diagram.  Specific  subroutines  for  plotting  the  information  are  not  shown  in  this  appendix. 


♦ 

*  THIS  PROGRAM  IS  THE  MAIN  PROGRAM  FOR  DATA  ACQUISITION 

*  OF  DATA  FROM  STORED  TAPES  OF  UNDERWATER  EXPLOSIONS 

*  VERSION  1 

*  LT  R.  A.  SHAFER.  USN 

*  NAVAL  POSTGRADUATE  SCHOOL 

*  28  APRIL  1987 

* 


*  SET  UP  VARIABLES 

*  MYBUFSZ  IS  THE  BUFFER  SIZE  FOR  TEMPORARY  DATA  STORAGE  OF  A/D  DATA 

*  FACTOR  IS  A  CONVERSION  FACTOR  FOR  BUFFER  MANAGMENT 

*  N  BUFFS  THE  NUMBER  OF  BUFFERS  FOR  DATA  COLLECTION 

*  NCLKS  THE  NUMBER  OF  CLOCKS  USED  FOR  A/D  CONTROL 

*  NEAR  FREQ  PARAMETER  FOR  CLOCK  CONTROL 

*  RDONLY  SET  A/D  CONVERSION  TO  READ  ONLY 

*  STARTLO  SET  CLOCK  WAVEFORM  IN  LOW  (ZERO  VOLTAGE)  POSITION 

*  TIMEOUT  TIME  LIMIT  FOR  DATA  ACQUISITION 

*  FREQ  FREQUENCY  OF  SAMPLING  CLOCK 

*  WIDTH  WIDTH  OF  SAMPLE  PULSE  (0.5  MICROSEC  MIN) 

INTEGER  MYBUFSZ 

INTEGER  FACTOR 

INTEGER  NBUFFS 

INTEGER  NCLKS,  NEARFREQ,  RDONLY 

INTEGER  RDWRT,  SLICE,  SQUARE 

INTEGER  STARTLO,  TIMEOUT 
REAL  FREQ,  WIDTH 

PARAMETER  (  MYBUFSZ  *  16384  ) 

PARAMETER  (  FACTOR  =  2) 

PARAMETER  (  NBUFFS  =  3  ) 

PARAMETER  (  NCLKS  =  1  ) 

PARAMETER  (  NEARFREQ  =  0  ) 

PARAMETER  (  RDONLY  =  1  ) 

PARAMETER  (  RDWRT  =  0  ) 

PARAMETER  (  SLICE  =  10  ) 
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PARAMETER  (  SQUARE  =  4  ) 

PARAMETER  (  STARTLO  =  0  ) 

PARAMETER  ( TIME0000  ) 

PARAMETER  (  WIDTH  =  0.0  ) 

SET  UP  INPUT  BUFFER  LIMITS 

ENSURE  BUFFER  IN  EVEN  WORD  BOUNDRY  IN  STACK 

INTEGER*2  DATABLK(98304) 

INTEGER  DUMMY 
EQUIVALENCE  (DAT  ABLK, DUMMY) 

DECLARE  SERVICE  ROUTINE  AS  EXTERNAL  PROGRAM 

EXTERNAL  ADSERVICE 

SET  UP  VARIABLES  FOR  A/D  CONTROL  AND  SETTINGS 

ADPATH  MARKER  FOR  PATH  TO  A/D  CONVERTER 

CLKPATH  MARKER  FOR  PATH  TO  CONTROL  CLOCK 

NCHAN  NUMBER  OFA  CHANNELS  TO  BE  SAMPLED 

FCHAN  FIRST  CHANNEL  TO  BE  SAMPLED 

SCHAN  SPACING  OF  CHANNELS  IF  MORE  THAN  ONE  SAMPLED 

GAIN  GAIN  SETTING  OF  A/D  CONVERTER 

INDEX  ARRAY  POINTER  FOR  BUFFER  MANAGEMENT 

POWER  NUMBER  TO  RAISE  2  TO  THE  POWER  OF  FOR  SAMPLE  BLOCK  SIZE 

TYME  ARRAY  FOR  TIME  VALUES 

RECLEN  TIME  RECORD  LENGTH  OF  DATA  SET 

RFREQ  ACTUAL  FREQUENCY  CLOCK  SET  TO 

R WIDTH  ACTUAL  WIDTH  OF  SAMPLE  PULSE 

SAMPLES  ARRAY  FOR  TIME  RECORD  DATA 

LOW  LOWER  LIMIT  FOR  DATA  STORAGE  ON  HARD  DISK 

HIGH  UPPER  LIMIT  FOR  DATA  STORAGE  ON  HARD  DISK 


INTEGER 

ADPATH 

INTEGER 

CLKPATH 

INTEGER 

STATWDS(2) 

INTEGER 

I 

INTEGER 

NCHAN 

INTEGER 

FCHAN 

INTEGER 

SCHAN 

INTEGER 

GAIN 

INTEGER 

INDEX 

INTEGER 

POWER 

REAL 

TYME(16384) 

REAL 

RECLEN 

REAL 

RFREQ 

REAL 

RWIDTH 

REAL 

SAMPLES!  16384) 

REAL 

LOW 

REAL 

HIGH 

REAL 

TIMSCL 

REAL 

SUM 

CHARACTER  ANSWER 


COMMON  SAMPLES,  ADPATH,  CLKPATH,  INDEX,  DATABLK,  NITEMS 


*  WRITE  HEADER  ON  SCREEN 

* 

WRITE  (6,9999) 

9999  FORMAT  (47H1  UNDERWATER  EXPLOSIONS  DATA  ACQUISITION  PROGRAM) 

* 

*  GET  INPUT  PARAMETERS  FROM  SCREEN 

♦ 

PRINT  *, "INPUT  NUMBER  OF  CHANNELS  TO  BE  SAMPLED" 

READ  *,  NCHAN 

PRINT  *, "INPUT  FIRST  CHANNEL  TO  BE  SAMPLED" 

READ  *,  FCHAN 

PRINT  *, "INPUT  SPACING  OF  CHANNELS  (  0  FOR  NONE  ) " 

READ  *,  SCHAN 

PRINT  INPUT  GAIN  FACTOR  (  0,1, 2, 3  =  *1,  *2,  *4,  *8  )'' 

READ  *,  GAIN 

* 

*  ENSURE  ALL  DACP  DEVICES  ARE  CLOSED 

CALL  MRCLOSALL 

* 

*  OPEN  A/D  CONVERTER  PATH 

*  SET  FIRST  CHANNEL,  SPACING  AND  GAIN  ON  A/D  CONVERTER 

* 

ADPATH  =  -l 

CALL  MROPEN(ADPATH,  "/dev/dacpO/adfO”,  RDONLY) 

CALL  MRADMOD(ADPATH,  0,  0) 

CALL  MRADINC(  ADPATH,  FCHAN NCHAN,  SCHAN  .GAIN) 

* 

*  IDENTIFY  BUFFERS  FOR  DATA  COLLECTION  AND  RELEASE  THEM 

*  TO  THE  EMPTY  BUFFER  QUEUE  OF  A/D  CONVERTER 

* 

DO  200  I  =  1.  81921,  MYBUFSZ 

CALL  MRBUFID  (ADPATH,  DATABLK(I),  MYBUFSZ*2) 

CALL  MRBUFREL  (ADPATH,  DATABLK(I)) 

200  CONTINUE 
* 

*  OPEN  CLOCK  PATH  FOR  A/D  CONVERTER  CONTROL  CLOCK 

*  AND  INTERNALLY  CONNECT  TO  THE  A/D  CONVERTER 

* 

CLKPATH  =  -1 

CALL  MROPEN  (CLKPATH,7dev/dacpO/efcUc2",  RDWRT) 

CALL  MR  WIRE  (ADPATH,  0) 

*  GET  CLOCK  SETTINGS  FROM  SCREEN 

* 

PRINT  *,  "INPUT  DESIRED  FREQUENCY  IN  HZ" 

READ  *,  FREQ 

* 

*  DISARM  CLOCK  TO  ENSURE  IT  IS  NOT  RUNNING 

*  AND  SET  UP  CLOCK  PARAMETERS 

CALL  MRCLKDIS  (  1,  CLKPATH  ) 
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CALL  MRCLK1GATED  (  CLKPATH,  NEARFRQ,  FREQ,  RFREQ,  1,  1.5,  R WIDTH, 
STARTLO,  5) 

* 

*  GET  NUMBER  OF  DATA  POINTS  TO  BE  RECORDED 

*  AND  TIME  SCALING  FACTOR 

PRINT  *, "INPUT  POWER  OF  TWO  FOR  NUMBER  OF  SAMPLES  (  MAX  =  14  )" 
READ  *,  POWER 
NITEMS  =  2“POWER 

PRINT  ‘."NUMBER  OF  SAMPLES  IS  "NITEMS 

RECLEN  =  NITEMS  /  RFREQ 

PRINT  *, "ENTER  TIME  SCALING  FACTOR" 

READ  ‘.TIMSCL 

PRINT  ‘."SIMULATED  SAMPLING  FREQUENCY  IS  ->  " ,RFREQ*TIMSCL 
PRINT  ‘."TIME  RECORD  LENGTH  IS  ".RECLEN/TIMSCL 
PRINT  ‘."START  TAPE  NOW" 

PRINT  ‘."THEN  ENTER  "y"  AND  HIT  RETURN  WHEN  TAPE  IS  READY" 

READ  *,  ANSWER 

* 

*  ARM  CLOCK  AND  WAIT  FOR  TRIGGER  TO  DROP  TO  ZERO 

* 

PRINT  ‘."ARMING  CLOCK” 

CALL  MRCLKARM  (  1.  CLKPATH  ) 

PRINT  *, "WAITING  FOR  TRIGGER  PULSE" 

* 

*  COLLECT  DATA  SET  AND  SEND  BUFFER  TO  A/D  SERVICE  ROUTINE 

CALL  MRXINQ  (  ADPATH,  NITEMS,  NITEMS,  ADSERVICE  ) 

CALL  MREVWT  (  ADPATH,  STATWDS,  30000  ) 

* 

*  COMPUTE  AVERAGE  OF  FIRST  1200  SAMPLES  TO  REMOVE  ANY  ZERO  OFFSET 

* 

PRINT  *,  "STOPPED  READING" 

DO  250  1=  1,  1200 
SUM=  SUM  +  SAMPLES(I) 

250  CONTINUE 

SUM  =  SUM/ 1200 

* 

*  FILL  TIME  ARRAY  AND  REMOVE  ZERO  OFFSET  FROM  RECORD 

* 

DO  300  I  =  1,  NITEMS 
TYME(I)=(I-  1)/(RFREQ*TIMSCL) 

SAMPLES(I)  =  SAMPLES  (I)  -  SUM 
300  CONTINUE 

9990  FORMAT  (IX,  G20.8,  10X  ,  G20.8) 

♦ 

*  CALL  PLOTTING  ROUTINE  FOR  DISPLAY  OF  DATA 

* 

CALL  PLOTl(  NITEMS,  SAMPLES,  TYME,  ”  Time  History  " ,  "Voltage"  ,  "Secon 

ds'V'timhis.g" ) 

* 

*  WRITE  DATA  TO  HARD  DISK 

* 

*  PRINT  *, "  INPUT  LIMITS  OF  HARDDISK  DATA  RECORD  " 

*  READ  *,  LOW.HIGH 

*  DO  350  I  =  1,  NITEMS 


*  IF  (  TYME(I)  .GE.  LOW  .ME(I)  .LE.  HIGH)  THEN 

*  WRITE  (7,9990)  TYME(I),  SAMPLES(I) 

*  END  IF 
*350  CONTINUE 

STOP 

END 

SUBROUTINE  ADSERVICE 
INTEGER  K 
INTEGER  INDEX 
INTEGER  POWER 
INTEGER  ADPATH 

INTEGER  NITEMS 
INTEGER*2  DATABLK(98304) 

INTEGER  GETINDEX 
PARAMETER  (GETINDEX  =  3) 

REAL  DUMMY 

REAL  SAMPLES(  16384) 

COMMON  SAMPLES,  ADPATH, CL KPATH,  INDEX,  DATABLK,  NITEMS 

IF  (.EQ.  0)  THEN 

NITEMS  =  16384 
END  IF 

*  STOP  CLOCK,  GET  BUFFER  FROM  DONE  BUFFER  QUEUE 

*  CORRECT  DATA  BY  SCALING  FACTOR  (20/2A12) 

*  FILL  SAMPLE  ARRAY  WITH  TIME  RECORD  AND  RELEASE  BUFFER  TO 

*  READY  BUFFER  QUEUE 

CALL  MRCLKDIS(  1,  CLKPATH  ) 

CALL  MRBU  GETINDEX,  INDEX) 

TEMP  =  I  /  204.8 
DO  300  K  =  1,  NITEMS 
DUMMY  =  FLOAT(DATABLK(K))  *  TEMP 
SAMPLES(K)  =  DUMMY 
300  CONTINUE 

CALL  MRBUFREL  (ADPATH,  DATABLK(INDEX)) 

RETURN 

END 
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APPENDIX  C 

IMPACT  TESTING  DATA  ACQUISITION  AND  REDUCTION  PROGRAM 

This  section  contains  the  program  listing  used  to  acquire  the  data  from  the  aluminum 
plate.  All  routines  were  written  by  the  author  except  where  noted.  For  further  explanation 
of  the  data  acquisition  specific  routines  the  reader  is  encouraged  to  refer  to  the 
MASSCOMP  documentation  "Data  Acquisition  User's  Manual".  Figure  47  is  the  program 
flow  diagram.  Sections  of  the  program  follow  the  general  outline  presented  in  the  flow 
diagram.  Specific  subroutines  for  plotting  the  information  are  not  shown  in  this  appendix. 


*  THIS  PROGRAM  IS  THE  MAIN  PROGRAM  FOR  DATA  ACQUISITION 

*  AND  REDUCTION  OF  IMPACT  TESTING  OF  AN  ALUMINUM  PLATE 

* 

*  R.  A.  SHAFER,  LT,  USN 

*  NAVAL  POSTGRADUATE  SCHOOL 

*  MONTEREY.  CA  93940 

*  27  MARCH  1987 

* 

* 


*  SET  UP  CONSTANTS  AND  VARIABLES 

* 


INTEGER 
INTEGER 
INTEGER 
INTEGER 
REAL 
REAL 
INTEGER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 
INTEGER*2 
INTEGER 
EQUIVALENCE 
EXTERNAL 
INTEGER 


MYBUFSZ,  NBUFFS 
NEARFREQ,  RDONLY 
RDWRT,  SQUARE 
STARTLO,  TIMEOUT 
FREQ,  WIDTH 
PI 

OUTFILE 

( NBUFFS  =  10  ) 

(  NEARFREQ  =  0  ) 

( RDONLY  =  1  ) 

( RDWRT  =  0 ) 

( SQUARE  =  4  ) 

(  STARTLO  =  0  ) 

(  TIMEOUT  =  30000  ) 

(  WIDTH  =  0.0 ) 
DATABLK(  100000) 
DUMMY 

(DATABLK.DUMMY) 
ADSERVICE 
ADPATH 
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INTEGER  CLKPATH(2) 

INTEGER  STATWDS(2) 

INTEGER  I 

INTEGER  NCHAN 

INTEGER  FCHAN 

INTEGER  SCHAN 

INTEGER  GAIN 

INTEGER  INDEX 

INTEGER 
REAL  TPI 

REAL  WEGHT(513) 

REAL  TYME(1024) 

REAL  RECLEN 

REAL  RFREQ 

REAL  R  WIDTH 

REAL  XTR(1024),  Y1TR(1024) 

REAL  XTI(1024) 

REAL  Y1TI(1024) 

REAL  XX(  1024) 

REAL  YY  1(1024) 

REAL  XFR(I024),XFI(1024) 

REAL  Y 1  FR(  1024),  Y 1  FI(  1 024) 

REAL  GXX(513) 

REAL  GYY(513) 

REAL  GXYRE(513) 

REAL  GXYIM(513) 

INTEGER  NOSAMP 

REAL  DELTAFRQ(1024) 

REAL  DFRE.Q 

REAL  HFR(513),HFI(513) 

REAL  HF(513),HMR(513) 

REAL  HMI(513),HM(513) 

REAL  AVG(5) 

CHARACTER  ANSWER 

COMMON  /BCR/  XTR,  YlTR,  ADPATH,  DATABLK,  MYBUFSZ,  NCHAN,  AVG 
COMMON  /CCSE/  TYME,  DELTAFRQ,  OUTFILE 
OUTFILE  =  8 

* 

*  OPEN  FILES  FOR  OUTPUT  OF  DATA 

♦ 

OPEN  (7  JTLE="/usr/shafer/pIate/dala"  ,STATUS="NEW") 

♦ 

*  CALCULATE  HANNING  WINDOW  WEIGHTING  FACTOR 

* 

TPI  =  8.0  *  ATAN(l.O) 

PI  =  4.0  *  ATAN(i.O) 

TEMP  =  TPI  /  1025.0 
SCL  =  SQRT(2.0/3.0) 

DO  401=  1,512 

WEGHT(I)=  SCL  *  (1.0-  COS  (TEMP*  FLO  AT(I))) 

40  CONTINUE 

* 

*  SET  UP  CONSTANTS  FOR  A/D  CONVERTER  SETUP  AND  BUFFERS 

* 

MYBUFSZ  =  5000 
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NCHAN  =  2 
FCHAN  =  16 
SCHAN  =  1 
GAIN  =  0 
ADPATH  =  -1 


*  OPEN  A/D  CONVERTER  PATH  AND  SET  CONVERTER  MODE 

* 

CALL  MROPEN( ADPATH,  7dev/dacp0/adfl".  RDONLY) 

CALL  MR  AD  INC(  ADPATH  FCHAN,  NCHAN,  SCHAN,  GAIN) 

CALL  MRADMOD(ADPATH,  0,  0) 

* 

*  OPEN  CLOCK  PATHD  PULSES  AND  SET  CLOCK  MODE 

* 

*  SET  UP  CLOCK  FOR  A/D  CONVERTER 

* 

CLKPATH(l)  =  -1 
CLKPATH(2)  =  -1 

CALL  MROPEN  (CLKPATH(l),7dev/dacp0/efclk0",  RDWRT) 

CALL  MROPEN  (CLKPATH(2),7dev/dacpO/efclk4",  RDWRT) 

FREQ  =  5000. 

PRINT  V  SETTING  UP  CLOCKS” 

CALL  MRCLK2  (CLKPATH(1),CLKPATH(2),0,0JFREQ,RFREQ,0.200000..RWIDTH,2,0) 

CALL  MRCLKTRIG  (ADPATH,  2.CLKPATH) 

print  *," Actual  samp  frequency  per  channel  --  ”,RFREQ 

NITEMS  =  60000 

PRINT  ‘."NUMBER  OF  SAMPLES  IS  "NITEMS 

* 

*  CALCULATE  TIME  RECORD  LENGTH  AND  FREQUENCY  RESOLUTION 

♦ 

RECLEN  =  1024/RFREQ 

PRINT  ‘."TIME  RECORD  LENGTH  IS  ".RECLEN 

DFREQ=  1 /RECLEN 

DO  60  I  =  1,  1024 

TYME(I)  =  I  /  RFREQ 

DELTAFRQ(I)  =  (I-1)*DFREQ 

60  CONTINUE 
* 

*  GET  NUMBER  OF  AVERAGES  TO  BE  COMPUTED  FROM  USER 

* 

PRINT  *, "INPUT  NUMBER  OF  AVERAGES  TO  BE  COMPUTED" 

READ  *  NOS  AMP 

* 

*  SET  UP  BUFFERS  AND  RELEASE  THEM  TO  A/D  PATH 
DO  200  I  =  1,  45001,  MYBUFSZ 

CALL  MRBUFID  (ADPATH,  DATABLK(I).  MYBUFSZ»2) 

CALL  MRBUFREL  (ADPATH,  DATABLK(I)) 

200  CONTINUE 

* 

*  START  DATA  ACQUISITION 

* 

DO  250  L  =  1  ,  NOSAMP 
65  PRINT  *,  L 
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PRINT  *,"A  WAITING  TRIGGER  ARM,  TYPE  ANY  CHARACTER  AND  RETURN 
WHEN  READY" 

READ  *,  ANSWER 
AVG(l)  =  0 
AVG(2)  =  0 

* 

*  TAKE  IN  DATA  RECORD 

* 

CALL  MRXINQ  (ADPATH,  MYBUFSZ,  NITEMS,  ADSERVICE)  MREATH.S,  30000  ) 
PRSTOPPED  READING” 

* 

*  DISPLAY  DATA  RECORD  FOR  CORRECTNESS 

*  IF  CORRECT  CONTINUE,  IF  NOT  DISCARD  DATA  SET  AND  TAKE  ANOTHER 

* 

CALL  PLOT1(1024,  XTR,TYME,"DATA  PREVffiW","VOLTAGE","TIME","pre.g") 

PRINT  *  ,  "Good  data  set?” 

READ  *,  ANSWER 
IF  (ANSWER  .EQ.  "n")  THEN 
GOTO  65 
END  IF 

* 

•  BEGIN  DATA  REDUCTION 

* 

DO  70  I  =  1.  1024 
XX(I)  =  XTR(I) 

YY1(I)  =  Y1TR(I) 

70  CONTINUE 

* 

*  WEIGHT  TIME  RECORD  WITH  HANNING  WINDOW  TO  ENSURE  PERIODICITY 

*  REQUIREMENTS  OF  FFT  ARE  MET.  XX.  BEING  AN  IMPACT  MEASUREMENT 

*  IS  PERIODIC  BY  DEFINITION  AND  REQUIRES  NO  WEIGHTING 

* 

DO  80  1=  1,512 

ITMP  =  1025  -  I 

YYl(I)  =  YYl(I)  *  WEGHT(D 

YYl(ITMP)  =  Y Y 1  (ITMP)  *  WEGHT(I) 

80  CONTINUE 

* 

*  CALCULATE  POWER  SPECTRA  INFORMATION 

*  FOLLOWING  PORTION  TAKEN  FROM  CCSE  PROGRAM 

* 

CAUL  FFT842(0, 1024.XX, YY 1  .OUTFTLE) 

GXX(l)  =  GXX(l)  +  4.0  *  XX(1)**2 
DO  90  K  =  2,  513 
J  =  1026  -  K 

GXX(K)  =  GXX(K)  +  (XX(KHXX(J))«2  +  (YY1(K)-YY1(J))**2 
GYY(K)  =  GYY(K)  +  (YY1(X)+YY1(J))**2  +  (XX(J)-XX(K))**2 
GXYRE(K)  =  GXYREIK)  +  XX(K)*YY1(J)  +  XX(J)*YY1(K) 

GXYIMOO  =GXYIM(K)  +XX(I)**2  +YY1(J)**2  -XX(K)**2  -YY1(K)**2 
‘Hi  CONTINUE 

GYYII)  =  GYY(I )  +  4  0*  YY I(  1  )**2 
GX  YRE(  1 )  =  GXYREI 1 )  ♦  2.0N  XX(  1  )* YY l<  1 » 

GXYIM(l)  =  0.0 


CALI  ZERO*  XT!  |n;-» . 

CALI.  ZEROtYin 

*  CALCULATE  FREQUENCY  RESPONSE  <>E  INPC7  xMx»(  ’!>!  ' 

*  SEPARATE  OF  THE  POWER  SPECTRA 

DO  100  I  =  1.  1024 
XX(I)  =  XTRih 
VYliI)  =  YITR.Ii 
100  COSTIVE 

*  WEIGHT  TIME  RECORD  WITH  HANNING  WINDOW  TOINs:  R(  E!  SO  •'  !<  !  ’ '• 

*  REQUIREMENTS  OF  FFT  ARE  MET  XX  BEING  AN  IMF’V  I  MI  W  Pi  Mf  N  ' 

*  IS  PERIODIC  BY  DEFINITION  AND  REQUIRES  NO  WEIGHTING, 

DO  1 10  I  =  1.  512 
ITMP  =  1025  •  I 
YYKI)  =  YYl(I)  *  WEGFTT(I) 

YYKITMP)  =  YY  liITMP)  *  WEGHT.Ii 
110  CONTINUE 

CALL  FFT842(0.1024.XX.XTI.OUTFILEi 
CALL  FFT842(0,1024.YYl.YITI.OLTFILEi 
TEMP  =  1 /FLOAT*  NOS  AMP) 

DO  1201  =  1.  513 
XFR(I)  =  XFR(I)  ♦  XX(I)*TEMP 
XFI(I)  =  XFI(I)  XTI(I)*TEMP 
Y1FR(I)  =  Y1FR(I)  +  YYKICTEMP 
YIFI(I)  =  Y1FKE)  +  YITKH'TEMP 
120  CONTINUE 

* 

*  TAKE  NEXT  DATA  SET 

♦ 

250  CONTINUE 

* 

*  WRITE  DATA  TO  HARD  DISK 

* 

DO  130  I  =  1.  513 

WRITE(7,9990)  GXX(I),GYY1I),GXYR£(I),GXYIM(I> 

130  CONTINUE 

9990  FORMAT(lX,G16  10.lx.G16  10.U.G16  10.UG16  lOi 

* 

*  WEIGHT  POWER  SPECTRA  INFORMATION  TO  CORRECT  FOR  V  MBE.R  <  >1  DM  \ 

SETS 

* 

CALL  COHERI 1 024 JCFREQNOS AMP.G  X X  ,G  Y  Y  G  X YRE  G  X  Y I M  I  u  I 

* 

*  CALCULATE  TRANSFER  FUNCTION  FROM  POWER  SPECTRA  INI  '  iRMA  1 1<  -N 

* 

DO  300  I»  1.512 
K  *  1025  -  I 

HFR(T)  =  GXYREH)/GXX(Ii 
HFR(K)  =  HFR(I) 

HFLI)  *  GXYIMIIiAjXXiIi 
HFRK)  *  HFKI) 

HFHn  10*  ALOG<SQRT(HFRlI,«2*Hnn.»*2  • 


HF<  K  i=  HFi  1 1 
;i<>  C  ONTINUE 

*  CALCl  LATE  THE  MOBILITY  TRANSFER  FUNCTION 

DO  350  I  =  2.  1024 

HMR(I>  =  HFRdi/DELTAFRQd) 

HMIdi  =  HFTd  i/DELTAFR  Q(I) 

HMi  It  =  10*  ALOG<SQRT(HMRd)**2+HMId)**2)) 

CONTINUE 

*  CALCULATE  THE  FREQUENCY  DOMAIN  OF  THE  TIME  SIGNALS 
DO  400  I  =  1 ,  5 1 3 

XFR(I)  =  SQRT(XFR(I)**2+XFI(I)**2) 

XFId)  =  10*  ALOG(XFR(I)) 

YlFRd)  =  SQRT(Y1FR(I)**2+Y1FI(I)**2) 

Y1FI(I)  =  10  *  ALOG(YlFRd)) 

400  CONTINUE 

*  DISPLAY  INFORMATION  ON  SCREEN  AND  PLOT  IF  DESIRED 

CALL  PLOTl(513.HFR,DELTAFRQ,"H(f)'7'Real","Freq  (Hz)","hfreq.g") 

CALL  PLOTl(513.HHX>ELTAFRQ/'H(D’7'Imag'\"Freq  (Hz)",”pfreq.g") 

CALL  PLOT ’^513, HFDELTAFRQ,"H(D”."db";’Freq  (Hz)”,"Hf.g") 

CALL  PLO'i  U513,HMR,DELTAFRQ, "Mobility", "Real'V'Freq  (Hz)",”Hft.g") 

CALL  PLOT1  (5 13.H  MI, DELTAFRQ, "Mobility", "Imag”,"Freq  (Hz)",  ’Hft.g ") 

CALL  PLOTl(513,XFR,DELTAFRQ,’X(0","Mag","Frcq  (Hz)”,"xf.g") 

CALL  PLOTl(513,XFT,DELTAFRQ,"X(f)","db","Freq  (Hz)’V'logxf.g") 

CALL  PLOTl(513,YlFR,DELTAFRQ,"Y(0","Mag  ’,"Freq  (Hz)","yf.g") 

CALL  PLOTl(513,YlFIJDELTAFRQ,"Y(0","db","Freq  (Hz)","logyf.g") 

STOP 

END 

*  SUBROUTINE  TO  TAKE  FILLED  BUFFERS  FROM  A/D  CONVERTER  AND 

*  CHECK  TO  SEE  IF  IMPACT  HAS  OCCURED.  IF  SO  REMOVE  DC  COMPONENT 

*  AND  TRANSFER  TO  ARRAY  FOR  FURTHER  REDUCTION  AND  RETURN  BUFFER 

*  TO  READY  BUFFER  QUEUE  AT  A/D  CONVERTER.  IF  NOT  RETURN  BUFFER 

*  TO  READY  BUFFER  QUEUE  AT  A/D  CONVERTER. 

SUBROUTINE  ADSERVICE 
INTEGER  K 

INTEGER  INDEX 

INTEGER  POWER 

INTEGER  ADPATH 

INTEGER  TEMP 

INTEGER  NCHAN 

INTEGER  MYBUFSZ 

INTEGER*2  DATABLK(  100000) 

INTEGER  GETINDEX 

PARAMETER  (GETINDEX  =  3) 

REAL  DUMMY 

REAL  AVG(5) 

REAL  DUMMY1 

REAL  THOLD 
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REAL  XTR(  1024) 

REAL  Y1TR(1024) 

COMMON  /BCR/  XTR,  Y1TR,  ADPATH,  DATABLK,  MYBUFSZ,  NCHAN,  AVG 
TEMP  =  ADPATH 

CALL  MRBUFGET  (1,  GETINDEX,  INDEX) 

» 

*  AVERAGE  DATA  SET  TO  DETERMINE  DC  OFFSET  OF  INFORMATION 

* 

IF  (AVG(l).EQ.  0.00)  THEN 

DO  250  K  =  INDEX  ,  INDEX  +  ( 99  *  NCHAN),  NCHAN 
AVG(l)  =  AVG(l)  +  FLOAT(DATABLK(K))/(2*20480) 

AVG(2)  =  AVG(2)  +  FLOAT(DATABLK(K+1))/(2*20480) 

250  CONTINUE 
END  IF 

* 

*  SET  THRESHOLD  TO  0.25V  ABOVFSET  OF  HAMMER 

* 

THOLD  =  AVG(l)  +  0.25 

DO  300  K  =  INDEX.  INDEX  +  MYBUFSZ,  NCHAN 
DUMMY  =  FLOAT(DATABLK(K))/(2*204.8) 

* 

*  IF  THRESHOLD  IS  EXCEEDED  GO  BACK  200  INDEXES  AND  RECORD 

*  INFORMATION  TO  ARRAYS  FOR  FURTHER  REDUCTION 

* 

IF  (  DUMMY  .GE.  THOLD)  THEN 
DO  350  I  =  -200 , 0,  NCHAN 
J  =  K  +  I 
L  =  (1/2)  +  101 

* 

*  ENSURE  BUFLAP  IS  CORRECT 

* 

IF  (J  LE.  0)  THEN 
J  =  J  +  100000 
END  IF 

DUMMY  =  FLOAT(DATABLK(J" 

DUMMY  1  =  FLOAT(DATABLK(J+l)) 

XTR(L)  =  (DUMMY/(2*204.8))-AVG(l) 

Y1TR(L)  =  (DUMMY  1/(2*204 ,8))-AVG(2) 

350  CONTINUE 
1  =  50 

* 

*  READ  IN  REMAINING  INFORMATION  INTO  DATA  SETS  FOR  REDUCTION 


DO  400  J  =  K,K  +  975*NCHAN,  NCHAN 
1  =  1+1 

DUMMY  =  FLOAT(DATABLK(J)) 
DUMMY  1  =  FLOAT(DATABLK(J+l)) 
XTR  (I)  =  (DUMMY/(2*204.8))-AVG(l) 
Y1TR(I)  =  (DUMMY  1/(2*204.8))- AVG(2) 
400  CONTINUE 


*  SET  INDEX  TO  EXIT  DATA  TRANSFER  LOOP 


K  =  INDEX  +  MYBUFSZ  +  1 
END  IF 
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non  ononn  o o n n n n n n n n n n n n o n n n n o o n 


300  CONTINUE 

* 

*  RELEASE  BUFFER  TO  READY  BUFFER  QUEUE  FOR  A/D  CONVERTER 

* 

ADPATH  =  TEMP 

CALL  MRBUFREL  (TEMP,  DATABLK(INDEX)) 

RETURN 
END 

SUBROUTINE  COHER(NNN,ISR.NDSJP,GXX.GYY,GXYRE,GXYIM,SFX.SFY» 

C 

C - 

MAIN  PROGRAM:  A  COHERENCE  AND  CROSS  SPECTRAL  ESTIMATION  PROGRAM 
AUTHORS:  G.  C.  CARTER,  J.  F,  FERRIE 

NAVAL  UNDERWATER  SYSTEMS  CENTER 
NEW  LONDON,  CONNECTICUT  06320 
MODIFIED  BY  R.A.  SHAFER.  LT  USN 

NAVAL  POSTGRADUATE  SCHOOL 
MONTERY  CA 
APRIL  1987 

TO  MEET  THE  REQUIREMENTS  FOR  ANALYSIS 
OF  REAL  TIME  DATA 

INPUT:  NNN  IS  THE  NUMBER  OF  DATA  POINTS  PER  SEGMENT 

4  <  NNN  <  1025 
ISR  IS  THE  SAMPLING  RATE 
NDSJP  IS  THE  NUMBER  OF  DISJOINT  SEGMENTS 
SFX  IS  THE  SCALE  FACTOR  FOR  THE  INPUT  DATA  STORED  IN 
THE  XX  ARRAY 

SFY  IS  THE  SCALE  FACTOR  FOR  THE  INPUT  DATA  STORED  IN 
THE  YY  ARRAY 


SPECIFICATION  AND  TYPE  STATEMENTS 
REAL  ISR 

REAL  XX(1024),  YY(1024) 

REAL  GXX(513),  GYY(513).  GXYRE(513),GXYIM(513) 

REAL  WEGHT(513),  PHI(513) 

REAL  LINE(50) 

REAL  DELTAFRQ(1024),TYME(1024) 

EQUIVALENCE  (WEGI(l)) 

COMMON  /CCSE/  TYME,  DELTAFRQ,  OUTFILE 

NNN  IS  THE  NUMBER  OF  DATA  POINTS  PER  SEGMENT 

ISR  IS  THE  SAMPLING  RATE 

NDSJP  IS  THE  NUMBER  OF  DISJOINT  SEGMENTS 

SFX  AND  SFY  ARE  SCALE  FACTORS  FOR  THE  INPUT  DATA 


NFFTS  =  NDSJP 
SMALL  =  .IE-20 


CALCULATE  CONSTANTS 


TPI  =  8.0*ATAN(1.0) 
DEG  =  360.0/TPI 


non  _  non  _  ^  non 


VARX  =  0.0 
VARY  =  0.0 
DT  =  1.0/  FLOAT(ISR) 

SF  =  SQRT(ABS(SFX*SFY)) 

NPFFT  =  NNN 
NNNPl  =  NNN  +  1 
NNND2  =  NNN/2 
NND21  =  NNND2  +  ! 

NP2  =  NPFFT  +  2 

ND2  =  NPFFT/2 

ND2P1  =  ND2  +  1 

DF  =  1.0/(DT*FLOAT(NPFFT>) 

FNYQ  =  FLOAT(ISR)/2.0 
CONST  =  0.25*DT/FLOAT(NNN) 

FLOW  =  0.0 

FHIGH  =  FNYQ 

ISTRT  =  IFIX(FLOW/DF)  +  1 

ISTOP  =  IFIX(FHIGH/DF)  ♦  1 

NORMALIZE  ESTIMATES 

FNSG  =  FLOAT(NFFTS) 

OFNSG  =  L0/FNSG 
TEMPI  =  CONST*OFNSG*SFX 
TEMP2  =  CONST*OFNSG*SFY 
TEMP4  =  CONST*OFNSG*SF 
TEMP3  =  2.0*TEMP4 
DO  90  K=1  J9D2P1 
GXX(K)  =  GXX(K)*TEMP1 
GYY(K)  =  GYYOO*TEMP2 
GXYRE(K)  =  GXYRE(K)*TEMP3 
GXYIM(K)  =  GXYIM(K)*TEMP4 
>0  CONTINUE 
VARX  =  0.0 
VARY  =  0.0 
DO  100  K=1,ND2P1 
VARX  =  VARX  +  GXX(K) 

VARY  =  VARY  +  GYY(K) 

00  CONTINUE 
VARX  =  VARX*DF*2.0/SFX 
VARY  =  VARY*DF*2.0/SFY 

CONVERT  GXX  TO  DB  AND  PLOT 

DO  110  1=1  JND2P1 
XX(I)  *  GXX(T) 

PHI(I)  *  10.0*  ALOG  10(  AMAX 1  (GXX(I), SMALL)) 

10  CONTINUE 

CALL  PLOTl(ND2Pl PHl.DELTAFRQ,  "Gxx  (QVdbYFREQ  (Hz)","loggxx.g  ’) 

CONVERT  GYY  TO  DB  AND  PLOT 

DO  140  I=IND2P1 
XX(I)  =GYY(I) 

PHKI)  *  1 0.0* AL OG  10( AMAX  1(GYY( I), SMALL)) 


nonno  n  r>. 


140  CONTINUE 

CALL  PLOTl(ND2Pl.  PHI,  DELTAFRQ  GyvifV.  dh  FREQ  '  H/i  ,  tw 
C 

C  COMPUTE  CROSS  SPECTRUM  AND  MAGNITUDE  SQUARED  COHERENT  I 
C 

IX)  2)0  K=I.ND2P1 

PHD  K)  =  GXYRE(Ki**2  +  GXYIM(K)**2 
X  X  ( K  i  =  SQRT(ABS<PHI(K)/<GXX(Ki*GYY(K)>)) 

2  '0  CONTINUE 

CAI  L  PI.OTIiND2PI  .XX. DEI  TAFRQ.  Coherence  .  ,  Freq  i  H/i  .  coher 

C 

C  TERMINATE  PROGRAM 

C 

RETURN 

END 

C - - 

C  SUBROUTINE:  ZERO 

THIS  SUBROUTINE  STORES  IN  A  FLOATING  POINT  ARRAY 


SUBROUTINE  ZERO( ARRAY.  NUMBR) 

INPUT:  ARRAY  =  AN  ARRAY  OF  FLOATING  POINT  VALUES  TO  BE 
ZERO  FILLED 

NUMBR  =  NUMBER  OF  ARRAY  VALUES 

DIMENSION  ARRAY(l) 

C 

DO  10  K=l. NUMBR 
ARRAY(K)  =  0.0 
10  CONTINUE 
RETURN 
END 
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Figure  47:  Impact  Testing  Data  Acquisition  and  Reduction  Flow  Diagram 
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APPENDIX  D 


STEADY  STATE  VIBRATION  DATA  ACQUISITION  AND  REDUCTION 

PROGRAM 


This  Appendix  list  the  computer  program  utilized  to  acquire  the  data  for  steady  state 
vibrations  analysis.  All  routines  were  written  by  the  author  except  where  noted.  For 
further  explanation  of  the  data  acquisition  specific  routines  the  reader  is  encouraged  to  refer 
to  the  MASSCOMP  documentation  "Data  Acquisition  User's  Manual".  Figure  48  is  the 
program  flow  diagram.  Sections  of  the  program  follow  the  general  outline  presented  in  the 
flow  diagram.  Specific  subroutines  for  plotting  the  information  are  not  shown  in  this 
Appendix. 


*  THIS  PROGRAM  IS  THE  MAIN  PROGRAM  FOR  DATA  ACQUISITION 

*  FOR  STEADY  STATE  ANALYSIS  OF  A  VIBRATING  BEAM 


*  07  MAY  1987 


INTEGER  MYBUFSZ,  NBUFFS 

INTEGER  NEARFREQ,  RDONLY 

INTEGER  RDWRT,  SQUARE 

INTEGER  STARTLO,  TIMEOUT 

REAL  FREQ,  WIDTH 

REAL  PI 

INTEGER  OUTFELE 

PARAMETER  (  NBUFFS  =  5  ) 
PARAMETER  (  NEARFREQ  =  0  ) 
PARAMETER  (  RDONLY  =  I ) 
PARAMETER  (  RDWRT  =  0  ) 
PARAMETER  (  SQUARE  =  4  ) 
PARAMETER  (  STARTLO  =  0  ) 
PARAMETER  (  TIMEOUT  *  30000  ) 
PARAMETER  (  WIDTH  =  0.0  ) 
INTEGER*2  DATABLK(40960) 

INTEGER  DUMMY 
EQUIVALENCE  (DATABLK, DUMMY) 
EXTERNAL  ADSERVICE 
INTEGER  ADPATH 

INTEGER  CLKPATH(2) 

INTEGER  STATWDS(2) 

INTEGER  I 
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INTEGER  NCHAN 

INTEGER  FCHAN 

INTEGER  SCHAN 

INTEGER  GAIN 

INTEGER  INDEX 

INTEGER  POWER 

REAL  TPI 

REAL  WEGHT(1025) 

REAL  TYME(2048) 

REAL  RECLEN 

REAL  RFREQ 

REAL  RWIDTH 

REAL  XTR(2048),  Y1TR(2048) 

REAL  Y1TI(2048) 

REAL  XX(2048) 

REAL  YY  1(2048) 

REAL  XFR(2048),XFI(2048) 

REAL  Y1FR(2048),Y1FT(2048) 

REAL  GXX(1025) 

REAL  GYY(1025) 

REAL  GXYRE(1025) 

REAL  GXYIM(1025) 

REAL  HWIN(2048) 

REAL  DELTAFRQ(2048) 

REAL  DFREQ 

REAL  MAG  1(1025) 

REAL  MAG2(1025) 

REAL  MAG3(1025) 

REAL  HFR(2048),HFI(2048) 

REAL  HR(2048),HI(2048) 

REAL  AVG(5) 

REAL  PI2,  TEMP 

INTEGER  AVGSAMP 

CHARACTER  ANSWER 

COMMON  /BCR/  XTR,  Y1TR,  ADPATH,  DATABLK,  MYBUFSZ,  NCHAN.  AVG 
COMMON  /CCSE/  TYME,  DELTAFRQ,  OUTFTLE 

* 

*  OPEN  DATA  FILES  FOR  DATA  STORAGE 

* 

OUTFTLE  =  8 

OPEN  (7  Xile="/usr/shafer/beam/siniuldala") 

* 

*  CALCULATE  WINDOWS  AND  CONSTANTS 

* 

TPI  =  8.0  *  ATAN(l.O) 

PI  =  4.0  *  ATAN(l.O) 

TEMP  =  TPI  /  2049.0 
SCL  =  SQRT(2.0/3.0) 

DO  401=  1,1024 

WEGHT(I)=  SCL  *  (1.0-  COS(TEMP*FLOAT(I))) 

40  CONTINUE 

PI2  =  2.0  *  ATAN(l.O) 

TEMP  =  PI2  /  512.0 
DO  46  I  =  2,  1024 
HWIN(I)  =  2.0 


v;  n'Tjrv’fti''  v* 


46  CONTINUE 

DO  47  I  =  1025,2048 
HWIN(I)  =  0.0 

47  CONTINUE 
HWIN(l)  =  1.0 

* 

*  SET  UP  A/D  CONVERTER  AND  BUFFERS  FOR  TEMPORARY  STORAGE  OF 

DATA 

* 

MYBUFSZ  =  4096 
NCHAN  =  2 
FCHAN  =  16 
SCHAN  =  1 
GAIN  =  0 
ADPATH  =  -1 

CALL  MROPEN( ADPATH,  "/dev/dacpO/adf  1 ",  RDONLY) 

CALL  MRADINC(ADPATH,  FCHAN,  NCHAN,  SCHAN,  GAIN) 

CALL  MR ADMOD( ADPATH,  0, 0) 

* 

*  INDENTIFY  BUFFERS  AND  RELEASE  THEM  TO  READY  BUFFER  QUEUE 

* 

DO  200  I  =  1,  36865,  MYBUFSZ 
CALL  MRBUFID  (ADPATH,  DATABLK(I),  MYBUFSZ*2) 

CALL  MRBUFREL  (ADPATH,  DATABLK(I)) 

200  CONTINUE 

* 

*  SET  UP  CLOCK  FOR  A/D  CONVERTER 
» 

CLKPATH(l)  =  -1 
CLKPATH(2)  =  -1 

CALL  MROPEN  (CLKPATH(l).7dev/dacp0/efclk0".  RDWRT) 

CALL  MROPEN  (CLKPATH(2),7dev/dacp0/efclk4",  RDWRT) 

FREQ  =  1000. 

PRINT  ‘."SETTING  UP  CLOCKS" 

CALL 

MRCLK2(CLKPATH(1),CLKPATH(2),0,0,FREQ,RFREQ,0,200000.,RWIDTH,2,0) 

CALL  MRCLKTRIG  (ADPATH,  2.CLKPATH) 

* 

*  COMPUTE  FREQUENCY  RESOLUTION  AND  FILL  FREQUENCY  ARRAY 

*  DISPLAY  RESOLUTION  ON  SCREEN 


RFREQ  =  RFREQ 

PRINT  Actual  sampling  frequency  per  channel  -  ".RFREQ 
PRINT  *, "INPUT  NUMBER  OF  AVERAGES  TO  BE  COMPUTED" 
READ  ‘.AVGSAMP 
NITEMS  »  4096 
RECLEN  =  2048/RFREQ 

PRINT  ‘."TIME  RECORD  LENGTH  IS  PER  CHANNEL  IS  "RECLEN 
DFREQ=  1 /RECLEN 

PRINT  FREQUENCY  RESOLUTION  IS  ".DFREQ 

DO  601=  1,2048 
TYME(I)  =  I  /  RFREQ 
DELTAFRQ(I)  =  (I-1)*DFREQ 
60  CONTINUE 
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•>>  ■  Wa*7i  ■>  »>  »vi  tfhiA  J 


MAIN  LOOP  FOR  DATA  ACQUSITION  AND  PRELIMINARY  REDUCTION 


DO  250  L  =  1  ,  AVGSAMP 

CALL  MRXINQ  (ADPATH,  MYBUFSZ,  NITEMS,  ADSERVICE) 

CALL  MREVWT  (ADPATH,  STATWDS,  0  ) 

PRINT  *,  AVG(3) 

DO  70  I  =  1,  2048 
XX(I)  =  XTR(I) 

YY  1(1)  =  Y1TR(I) 

70  CONTINUE 

* 

*  WEIGHT  TIME  RECORDS  WITH  HANNING  WINDOW 

* 

DO  80  1=  1,  1024 

ITMP  =  2049  - 1 

XX(I)  =  XX(I)  *  WEGHT(I) 

YY1(I)  =  YY1(I)  *  WEGHT(I) 

XX(ITMP)  =  XX(ITMP)  *  WEGHT(I) 

YYl(ITMP)  =  YY1(ITMP)  *  WEGHT(I) 

80  CONTINUE 

* 

*  COMPUTE  FORWARD  FFT 

* 

CALL  FFT842(0.  2048,  XX,  YYl.OUTFILE) 

* 

*  COMPUTE  SPECTRA 

* 

GXX(l)  =  GXX(I)  +  4.0*XX(1)**2 
DO  90  K=2,  1025 
J  =  2050  -  K 

GXX(K)  =  GXX(K)  +  (XX(KHXX(J))**2  +  (YY1(K)-YY1(J))**2 
GYY(K)  =  GYY(K)  +  (YY1(K)+YY1(J))**2  +  (XX(J)-XX(K))**2 
GXYRE(K)  =  GXYRE(K)  +  XX(K)*YY1(J)  +  XX(J)*YY1(K) 

GXYIM(K)  =  GXYIM(K)  +XX(J)**2  +YY1(J)**2  -XX(K)**2  -YY1(K)**2 
90  CONTINUE 

GYY(l)  =  GYY(l)  +  4.0*YY1(1)**2 
GXYRE(l)  =  GXYRE(l)  +  2.0*(XX(1)*YYI(1)) 

GXYIM(l)  =  0.0 


250  CONTINUE 
♦ 

*  DATA  COLLECTION  COMPLETE,  COMPLETE  REDUCTION  AND  DISPLAY 

RESULTS 

* 

CALL  COHER  1(2048 JIFREQ^ VGS AMP, GXX,GYY,GXYRE,GXYIM,  1.0, 1.0) 

* 

*  CALCULATE  FREQUENCY  RESPONSE  (MOBILITY) 

* 

DO  600  I  =  2,  1024 
K  =  2049  - 1 

HFR(I)  =  (GXYRE(I)/GXX(I))/DELTAFRQ(I) 

HFR(K)  =  HFR(I) 


HFI(I)  =  (GXYIM(I)/GXX(0)/DELTAFRQO) 

HFI(K)  =  HFI(I) 

600  CONTINUE 

HFR(l)  =GXYRE(1)/GXX(I) 

HFR(2048)  =  HFR(l) 

HF1(1)  =  GXYIM(1)/GXX(1) 

HFI(2048)  =  HFI(l) 

* 

*  RECORD  POWER  SPECTRA  INFORMATION  TO  DISK 

* 

DO  700  I  =  1,  1025 

WRITE  (7.9990)  GXX(I),GYYt  I),GX  YREiI  ),GXY  IMlU 
700  CONTINUE 

9990  FORMAT  (IX,  G16.8.G16.8.G16.8,GI6.8> 

* 

*  PLOT  RESULTS 
» 

CALL  PLOT1(1024.HFT?.DELTAFRQ,  H(0".  Real".  Freq  (H/i  .  Hfr  n  > 

CALL  PLOTl(  1024, HF1.DELT AFRQ," H(0",”Imag".  Freq  <  Hz).  HI) .n  i 
STOP 
END 

SUBROUTINE  COHERKNNN.ISRN’DSJP.GXX.GYY.GXYRE  GXYIM.SFX.SFY 
C 

C - 

C  MAIN  PROGRAM:  A  COHERENCE  AND  CROSS  SPECTRAL  ESTIMATION  PROGRAM 
C  AUTHORS:  G.  C.  CARTER,  J.  F  FERRIE 

C  NAVAL  UNDERWATER  SYSTEMS  CENTER 

C  NEW  LONDON,  CONNECTICUT  06320 

C  MODIFIED  BY  R  A.  SHAFER.  LT  USN 
C  NAVAL  POSTGRADUATE  SCHOOL 

C  MONTERY  CA 

C  APRIL  1987 

C  TO  MEET  THE  REQUIREMENTS  FOR  ANALYSIS 

C  OF  REAL  TIME  DATA 

C  INPUT:  NNN  IS  THE  NUMBER  OF  DATA  POINTS  PER  SEGMENT 

C  4  <  NNN  <  1025 

C  ISR  IS  THE  SAMPLING  RATE 

('  NDSJP  IS  THE  NUMBER  OF  DISJOINT  SEGMENTS 

C  SFX  IS  THE  SCALE  FACTOR  FOR  THE  INPUT  DATA  STORED  IN 

C  THE  XX  ARRAY 

C  SFY  IS  THE  SCALE  FACTOR  FOR  THE  INPUT  DATA  STORED  IN 

C  THE  YY  ARRAY 

C - 

C 

C  SPECIFICATION  AND  TYPE  STATEMENTS 

C 

REAL  ISR 

REAL  XX(2048),  YY(2048) 

REAL  GXX0025).  GYY(1025),  GXYRE(1025),  GXYIM(1025) 

REAL  WEGHT(1025),  PHI(1025) 

REAL  LINE(50) 

REAL  DELTAFRQ(2048),TYME(2048) 

EQUIVALENCE  (WEGHT(l)RHI(l)) 

COMMON  /CCSE/  TYME,  DELTAFRQ,  OUTFILE 


non 


c 


C  NNN  IS  THE  NUMBER  OF  DATA  POINTS  PER  SEGMENT 
C  ISR  IS  THE  SAMPLING  RATE  ' 


C  NDS/P  IS  THE  NUMBER  OF  DISJOINT  SEGMENTS 
C  SFX  AND  SFY  ARE  SCALE  FACTORS  FOR  THE  INPUT  DATA 


VFFTS  =  NDSJP 
SMALL  =  IE-20 

C 

C  CALCULATE  CONSTANTS 

C 

TP  I  =  X.0*ATAN<  1.0) 

DEG  =  WI.O/TPI 
VARX =00 
VARY  =  0.0 
DT  =  1.0/  FLOAT(ISR) 

SF  =  SQRTi  ABS(SFX*SFYi) 
N'PFFT  =  NNN 
NNNPI  =  NNN  +  I 
NNND2  =  NNN/2 
NND21  =  NNND2  +  1 
NP2  =  N'PFFT  +  2 
ND2  =  NPFFT/2 
ND2P1  =  ND2  +  I 
DF  =  1  0/(  DT*  FLO  AT  fNPFFT) ) 
FN'YQ  =  FLOAT(ISR)/2.0 
CONST  =  0.25*DT/FLOAT(NNN) 
FLOW  =  0  0 
FHIGH  =  FNYQ 
ISTRT  =  [R  X(  FLOW/DR  +  1 
ISTOP  =  IFIX(FH1GH/DF)  ♦  1 


NORMALIZE  ESTIMATES 


FNSG  =  FLOAT(NFETS) 

OFNSG  =  1.0/FNSG 
TEMPI  =  CONST*OFNSG*SFX 
TEMP2  =  CONST*OFNSG*SFY 
TEMP4  =  CONST*OFNSG*SF 
TEMP3  =  2.0*TEMP4 
DO  90  K=1.ND2P1 
GXX(K)  =  GXX(K)*TEMP1 
GYY(K)  =  GYY(K)*TEMP2 
GXYRE(K)  =  GXYRE(K)*TEMP3 
GXYIM(K)  =  GXYIM(K)*TEMP4 
90  CONTINUE 
VARX  =  0.0 
VARY  =  0.0 
DO  100  K=1  J9D2P1 
VARX  =  VARX  +  GXX(K) 

VARY  =  VARY  +  GYY(K) 

100  CONTINUE 

VARX  =  VARX*DF*2.0/SFX 
VARY  =  VARY*DF*2.0/SFY 


o  o  o  tJ  non  non 


i  i  uN'.t-K’  i,\\  Tu:  H  V  V> 

:  .< )  1 1"  ui  \d:p. 

xxn  =  < , \ x  : 

f’Hii  i  ■  =  :<•  - 1*  i  k  . :  ..  \m  x  \  <  xx  •  v  v . 

i  |)ST!V  f 

«  xi  i  pi  I  >t  ■  nd;pi  pm;  :  -t ;  -  *r»v  . .  ■  - 

( 

-  r  nsvf  p'  i .%  >  •< ,  ;>h  w  >  !<•  • 

;  *  »  :  -ii 1  :  -  s ;  > '  i  ■  ■ 

xx  ;  =<->'!  ; 

PHI.  I •  -  pi  i'*-v;  m  . :■  .  \m  \ x  ' 

.  (  ( >V7  IM  F 

<  \i  i  ;m  i  )T :  % ; >.'p'  ph:  :  .f :  •  - 

t 

(  <  1  (Ml*!  TT  ■\V):>!S!>!  \>  \  .’U\  , 

( 

'AMMO  =m. 

i*HJ.  I  •  -  "  11 

do  ;u>  k=;.nd:pi 

XXK  =  f.XVRF  K  ■ 

IF  •  XX K.  I >xi  I'f  !■*«> 
ro  IF  ■  fiXVlMi K  "  1'#)  IHI)  p*. 
ixn  XXK  =  0 

I  v<>  PHI.  Ki  =  DEG*ATAN2.GXYIM.  K  XXK 
;ix)  CONTINUE 
r 

C  PLOT  PHASE  FROM  PHLIM  TO  PHUM 

C 

PHLIM  =  1800.0 
[XT  210  K=2.ND2PI 
X  =  PHI(K)  -  PHI.  K  1 1 

PHKK)  =  PHI(K)  -  SIGNO60 ..X  i*  AIN'T  0  ‘v  ABS  X  ,  '**< 

IF  (PHI(K).GT  PFILINf)  PHI(K)  -  PHI.K.  PHLIM 
IF(PHKK).LT.HPHLIMt)  PHIfKi  =  PHIiKi  *  PHLIM 
210  CONTINUE 

COMPUTE  CROSS  SPECTRUM  AND  MAGNITUDE  SQUARED  CoHFKFVf. 
DO  230  K=1,ND2P1 

PHI(K)  =  GXYRE(K)**2  +  GXYIM(K)**2 
XX(K)  =  SQRT(ABS(PHI(K)/(GXX(K>*GYYiK)») 

30  CONTINUE 

CALL  PLOTl(ND2Pl  .XX, DELTAFRQ, "Coherence", " ",  Freq  ih/i  ,  coher  g  . 
CONVERT  GXY  TO  DB  AND  PLOT 
DO  250  1=1  JND2P1 

PHI(I)  =  5.0*ALOG10(AMAX1(PHI(I).SMALL)) 

50  CONTINUE 

TERMINATE  PROGRAM 

RETURN 
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'  1  AIM*  iN  H<  W  \  ’W  ■  -  V  •  «!•.*• 

■ !  n u*  «  ’>< 

V  .  >  ,  ■  -  \  ’.  *  M  *  • 

Vvl,.!  =  " 

V*  ( ii  2  ■  *  11 


GET  DONE  HI  EELK  FKuM  D«>\I  HI  i  I  IK  'jl  I  I  I  <  >\  VI)  '  " 

CALL  MRBLEGET  .  I  GETINDEX  IN  LEX  ■ 

AVERAGE  ENTIRE  ARRAY 

DO  250  I  =  INDEX.  INDEX  ♦  MYBL'FSZ  ,  2 

AVOID  =  AVG(  1 )  +  FLOAT ( DAT ABLK< I ) »/( 2  •  4 1 94 30.4 ; 
AVG(2)  =  AVG(2)  ♦  FLOAT<DATABLK(M »/(2.  *  419430  4) 
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